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ABSa-RACT 


Stress intensity factors in opening and sliding mode 
nave been computed, with a finite element scheme using 
special crack tip element of large size and consequently a 
coarse mesh for thin orthotropic rectangular plates. 
Theoretically exact stress function has been used in the 
special crack tip element with a good number of terms in 
the series. The compatibility of displacement across the 
special element boundary has been satisfied in a least square 
sense . 

Results have been verified in isotropic case using the 
same scheme. The maximum error with the theoretical value is 
only 4%. Energy release rate has been computed successfully 
in a case v^here direct computation of stress intensity 
factors are not possible. Evaluation of stress intensity 
factors have been done in primarily three different cases. 

They are a thin plate with a (1) central crack, (2j symmetric 
notches on both sides and (3) a notch in one side, where 
the loading is uniform tensile stress. Besides this case (3) 
has been dealt with under uniform shear stress and also a 
combination of uniform tensile and shear stress. 

A really good convergence of the computed values of 
stress intensity factors with the increase in number of degrees 
of freedom in the special crack tip element has been observed 
in all the cases , The linearness of the relation of applied 
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stress and computed stress intensity factor values 
has been verified. Dependence of computed values of 
stress intensity factors on poisson’s ratio has been 
studied in all three cases -which unlike the isotropic 
case is unique in an orthotropic elastic field. 
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Elemental displacement vector 
Displacement vector for ith node 
Unknov;n degrees of freedom 
Stress intensity factor in mode I . 

Stress intensity factor in mode II 
• Stress intensity factor in mode III 
Strain shape function matrix 
Modulus of elasticity 

Percentage error with respect to Kj & values 
corresponding to N =16. 

Modulus of elasticity in principal material 
direction 1 . 

Modulus of elasticity in principal material 
direction 2. 

Nodal force vector 

Energy release rate 

Modulus of shear rigidity in 1 ~2 direction 
Global^ stiffness matrix in isotropic case 
Global stiffness matrix in orthotropic case 
Element stifness matrix in isotropic case 
Element stifness matrix in orthOtrppic case 
Transformation matrix 

Number of degrees of freedom excluding rigid 
body translation and rotation u j ^ ^ 
special crack tip element. 

Radial distance from the origine 

Radius of the special crack tip element 
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Modified elasticity matrix 
Geometric transformation matrix 
Displacement in X -directi on 

Nodal displacement of the ith node in X-direction 
Displacement in Y-direction 

Nodal displacement of the ith node in Y-direction 

Rotational displacement 

Strain energy 

Complex argument 

Complex argument 

Real part of the root of the characteristic 
equation 

Imaginary part of the root of the characteristic 
equation. 

Stress function 

The vector containing unkno-wn degrees of freedom 

^4n* ^o’ ^o’ 

Total potential of a structure 

Root of the characteristic equation 

Poisson's ratio 

Poisson's ratio in principal material direction 
1-2 

Normal stress 
Stress in X-direction 
Stress in Y-direction 

Applied normal s'tress 30 degree inclined to 
Y™di recti on 


- 


Shear stress in X-Y direction 


Note ; Other symhols used in this -work have been clarified 
in places where-ever they are used. 
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CHAPTER - 1 


INTRODUCTION 

The rapid cultivation .and development of fracture mechanics 
has been taking place cwing to the fact that crack is the 
root cause of structural failure in many cases, 'A majority 
of large scale failures in such diverse structures as storage 
tanks, pressure vessels, pipelines, bridges, turbine generator 
rotors, ships, aircrafts, rocket motors etc., have been 
invariably traced to the presence of a crack originating 
from a site of stress concentration. The famous accidents 
in the year of 1954 involving Comet” jet aircrafts on schedule 
flights were traced to fatigue cracks at the corners of the 
windows, due to stress concentration in the fuselage. 

Accounts can be given on many other structural disasters 
whose root cause lie in development of cracks. 

1 . 1 Fracture mechanics and stress intensity factor 

The stress and strain patterns change abruptly, when 
a crack is introduced in a solid body, in the regions near 
the crack tip. This phenomenon gradually diminishes as one 
moves further and further away from the crack tip. In such 
a cracked body, when external loads are applied, the stresses 
increase rapidly in the vicinity of the crack tip and become 
unbounded when one reaches the crack tip. Linear fracture 
mechaniesis based on the premise that the stress-state in the 
vicinity of the crack tip can be characterized by stress 
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intensity factors. There are distinctly three different 
modes of fracture and they are ; (1) opening mode, (2) sliding 
mode and (3) tearing mode as illustrated in Fig. 1 (a) . A 
general fracture can he any possible combination of these 
different modes of fracture. Conventionally, the stress 
intensity factors in these three modes are represented as 
Kj, Kjj and K^j-j respectively. 

The relations of stress and displacement components 

■with the stress intensity factor may be illustrated as follows. 

Consider a simplified case of a crack in a uniform tensile 

field, see Fig. 1(b). Here rj- . and in the vicinity 

- y j4y 

of crack tip can be expressed by the following expressions. 

COS (©/2) (1 ■“ Sin (0/2). sin (3 ©/3))(1.1) 

^ (2%r)^ 

. COS (0/2) (1 + sin (0/2). sin (3 0/2)) (1.2) 

y (27ir)2 


^ sin (©/2) ® cos (©/2) « cos (,3 ©/2)** 

(27tr)* 

The two displacement components u & v along x and y axes 
can be written by the following expressions. 

K-r 1 

u = (r/2 (2k-1) cos^e/2) - cos (3 @/2) (1 .A) 

Kt 1 

^ (2k+1) sin (0/2) » sin (3 0/2) (1.5) 

where k= (3'-4iJ }for plane strain and k=:(3"l) )/(1 + V ) for 
plane stress. 
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K i.e. the stress intensity factor in Mode I fracture 
I 

is only dependent on the geometry of the cracked body and 

magnitude and configuration of the applied loading. There 

is another important parameter used in fracture mechanics 

2 

called the energy release rate, denoted as 'G' . When a orack 

is extended hy a small length of 'da’ there is some work done 

by the external loads which is equal to the amount of energy 

released'! ' from the body. Therefore G = . To predict 

whether a crack will propagate in a body under a set of external 

load we need t^ know a critical value of G i.e. critical 

energy release rate 'G„' or a critical value of stress 

intensity factor which is otherwise called 'fracture toughness' 

and noted as K. . K. and G* are dependent on material 
c c c 

properties. V/hen the stress intensity factor or energy 

release rate, which are only dependent on loading configuration 

and geometry of the cracked body, surpasses the critical 

values K and G. respectively a stable or unstable crack 
c c 

propagation will take place. If the crack propagation is 

unstable it leads to failure inevitably. Therefore, when a 

cracked structural member is under service condition one 

should be careful about either the energy release rate or 

the stress intensity factor does not exceed the respective 

critical value , Elnergy release rate and stress intensity 

2 3 

■factors are related by the following relations. ’ 
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For Mode I fracture in plain stress ; G 
For Mode I fracture in plain strain i G 

For Mode II fracture in plain stress; G 
For Mode II fracture in plain strain; G 


I " 

T /-2 

^"I 

E 


(1 .6) 

I = 


E 

(1.7) 

II"^ 

^11 
E ■ 


(1 .8) 

II"" 

!L 

(1-1)^) 

E 

(1.9) 


1 .2 E_valuation of S tress_ i_nt^ _^ ty_J^_ac^^ ^ 

1 .2.1 F V aluation M e th ods 

Analytical attempts for elastic stress analysis of cracked 
todies ■were first made as early as in the year of 1913. It 
has de'veloped considerably ' in the last two decades. Numerous 
contributions have been made in the development of analytical 
and numerical techniques. As it stands today, the complex 
variable techniques can handle a variety of crack shapes 
while finite element methods have attained an impressive 
position for application to problems with complicated geometry. 
Recently people have started \ising experimental techniques 
for the, determination of stress intensity factors. Experimental 
methods are intended as alternatives to analytical techniques 
for practical complex configurations. But as a matter of 
fact the rapid progress of two-dimensional analysis with 
complex variable and finite element methods are at every 
stage pushing the ex'perimental methods into unattractive 
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positions. However, even today three dimensional analyses 
are analyticallly difficult and computationally so unwieldy 
and expensive that for three dimensional problems, expertnental 
methods are still preferred to analytical and numerical 
techniques. Amongst the experimental methods both two and 
three dimensional photoelastic techniques are widely used 
and their choice is based on their simplicity. 

1.2.2 ela st ic_ _f ieltl 

With exception to some of the recent advancements, specially 
in aerospace engineering, the conventional structural materials 
are all metals or its alloys. As a matter of consequence, 
almost all the developments in fracture mechanics have dealt 
vjith the isotropic elastic fields. It is worth mentioning here 
that till early seventies there has not been any considerable 
development made in anisotropic field so far as the linear 
fracture mechanics is concerned. It is only since mid- 
seventy that people have started realising the hidden 
potential of fibre composite materials specially in the fields 
where weight reduction is one of the most important criteria 
of structural design. Currently almost every aerospace 
company is developing different products made of fibre composites 
and trying to successfully use them in aircrafts to substitute 
c .nventional structural materials. It is presently in the 
experimental stage and companies like 'General Dynamics* and 
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Me Donnell Douglas are using fibre oomposites in tiorizontal 
and vertioal stabilisers , for some portions in the fuselage 
as \Mell as engines . Currently vjith various conventional 
metal alloys j thrust to -weight ratio of 5 to 1 has been 
achieved. Reinforced plastics may lead to thrust to weight 
ratio as high as 16 to 1 . Ultimately with advanced graphite 
fibre composite thrust to weight ratio as high as 40 to 1 
apioear possible However , the road to this goal can be 
perilous as the failure behaviour of fibre -composites are not 
as well known as that of metals Or metal alloys. Inspite of 
the past developments of fracture mechanics in isotropic 
field; the study of cracks in orthotropic or other -anisotropic 
materials are of immense current importance to cater to the 
practical necessities. Today lot of work is being carried 
out all over the world to study the different failure 
phenomena in orthotropic or other anisotropic elastic fields 
to develop a profound theoretical backgreund for failures in 
fibre composites. 

1 .3 Literatiye Su r^e_y 

In early days of development of finite element methods, 

i 

the crack problems have been dealt with conventional triangular 
and other higher order elements. It is known from the 
continuum analysis that a singularity in stresses of the 
order of r ^ exists at the tip of the crack. For these problems. 
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very soon, the inability of conventional elements to represent 
these situations of large stress concentration or to represent 
the regions of large stress gradients is fast realized. Using 
conventional elements the stress distribution derived are 
in substantial error in the immediate vicinity of the crack. 

The stress intensity factors v;ere determined only by extra- 
polating the s cresses a-way from the crack tip with the known 
forms of near-crack tip stress field. Considering the 
rectangular plate with the central crack under pure tension 

1 

the near tip stress -field is known of the form, the equations 

( 1 . 1 ), (1 . 2 ) & (1 . 3 ) . 

The stresses derived from the conventional finite element 
analysis slightly away from crack tip. were fitted with this 
distribution and was determined. However, this method has 
not been found to yield accurate estimate of stress intensity 
factor. Then the realisation came that a different method of 
introducing crack tip singularity is necessary into the 
finite element methods. This lead to the development of methods 
of analysis wherein the region close to the crack tip is 
treated with special singular elements and the regions away 
from the crack tip are ddalt with conventional finite elements. 
One of the significant development is the hybrid method where 
the displacement description in singular element is drawn from 
the continuum solutions, either full solution or sometime only 
the relevant part of the solutions representing the crack tip 
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elastic singularity. Extensive literature has appeared in the 

development of hybrid methods and significant part of it is 

4-1 3 

listed in the ref. . Besides hybrid method, the develop- 

ment of isoparametric elements leading to singularity of the 
order r”^ at the crack tip has been another major breakthrough 
in this areac With the development of isoparametric elements 
it is felt thao singular crack tip elements of hybrid type are 
not necessary. 


4 

Some special crack tip elements have been used by Byskov , 
5 -- 

and Tracey--^ wliich represent (1/r)^ stress singularity for 
elastic analysis. But the element displacement functions, 
hov/ever, did not satisfy inter-element compatibility criteria. 
Pian and Tong have developed stress hybrid model using the 


modified principle of minijnura complementary energy for which 


the functional is %. 


^ f S . . , 1 0^ . . O' 


ijkl ^ Lj ^ kl dv -/ 

6V, 


Tj_Vj^ ds+ 


/ii 


ds 


O .11 ) 


^ij = stress tensor 

^i^kl ~ compliance tensor 

Vm = volume of the element 

<5v^ “ boundary of the elem.ent 

S = part of on which tractions are specified 

Vj^ boundary displacements 

= boundary tractions . 
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The complimentary displacement hybrid model had been used 

7 8 

by Atluri' ’ using the functional : 


m 


®ia ^kl-Fi “ 1 )<SV -H / (V-u^) ds - ^T.V^ 


ds 


V, 


m 


V. 


m 


(1 . 12 ) 


■where 

^ijkl stiffness tensor 

Gfj = i ( u^/ X 0 + uo/ xi) (1.13) 


^ r 


.-ij 


interior displacements 


For singular elements, the displacements are assumed as 

^ I 


= (u„ ) + (U ) \ ^ s 

^ K . S s -J , 

J 

= (u^)’ ■iP^^(U3) (K3j 


^1.14) 


in ■which (Up ) is simple polynomials and (U„ ) , are known as 
displacement functions for plane problems with crack, obtained 
from crack tip s-tress fields. The inter-element boundary 
displacemient is assumed in terms of nodal displacements (q^^) 

i.e . 

(1.15) 

where is the transformation matrix. The inter element 

compatibility is thus assumed. On the element boundaries 
radiating from cracm tip (r)'=’ type displacement behaviour 
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is built in. For all elements around orack tip i.e. the 

stress intensity factor is same, and final equations are 

obtained in terms of nodal displacements q and the stress 

intensity factors . 

s j 

There have been a series of attempts around the same time 
in the development of hybrid elements. A hybrid method of 
analysis was developed at Indian Institute of Science by 
Rao, Ra^u and Krishna Murthy^ for general problems of stress 
concentration. It is based on the development of large 
primary elements in the regions of stress concentration and 
stiffness of the primary elements are determined from an 
assumed displacement pattern provided by the continuum 
solutions. The stiffness matrix is finally derived from mini- 
mization of potential energy. The outside region was filled up 
by conventional elements. Alternative hybrid formulations 
around the same time has been done by Morley, Pian°. Another 
development is by Benzley who undertook the development of a 
generalised quadrilateral finite element that includes a 
singular point at a corner. Interelement displacement continuity 
i.e. compatibility has been maintained so that convergence 
of the finite element is preserved. A global -local concept 
of finite element formulation is utilized to formulate the 
the complete set of stiffness relationships. 
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Tha superposition approacli to finite element fracture 

1 1 

analysis has been successfully applied by Yamamoto and 
similar idea has also been used by Morley. The super imposition 
approach to finite element fracture analysis attempts to 
determine stress intensity factors through a linear combination 
of olassical singularity solutions and a coarse finite -element 
grid. The concept involves first determining the classical 
solution for crack in a infinite elastic body as close to the 
problem of interest as possible. Now, the finite element method, 
because of its ability to model elastio bodies charaoterized 
by complicated displacement and force boundary conditions as 
well as oomplex geometries is called upon to provide the 
second solution. The results from the finite element solution 
and classical solution are superimposed to determine the final 
solution. 

1 2 

R.D.. Hens hell &. K.G. Shaw have made a remarkable 
contribution in this field by devising an isoparametric 
singular element. A standard eight noded isoparametric element 
in an X-Y space is transformed to a square in the'f-'^l space 
with vertices at ( +1 > +'l ) • The midside nodes are shifted 
from their original position towards the crack tip at quarter 
point to exhibit a singularity of the order of r”^. The 
considered interpolation along one edge say edge 1-2 (see E'ig.1.(c 
with nodes at -1 , 0 & +1. in the parameter plane. Then r is 

X - X. 

defined as r — 

X2- x^ 


5 along the edge 1-2 r varies from 
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0 to 1 5 and the corner nodes are at r=0 & r=1 and midside 

node is at r=p. The interpolation functions are assumed as 

2 

r = a^ -r a^ t -5 

2 

and u = + 1 ) 2 ^ + b^'f 

Then the relation bet-ween r and is obtained as : 

2 ^ 

;i jL A l"8 r -!■ 1 6 p^) 

"^r “ ^ (1'"8 p 16 -!- 8(1 “-Z pj r) ^ 

tends to infinity when r = ° This singularity 

occurs at r=0 i.e. node 1, when p=-i.i.e. when the mid side 

node is shifted to the quarter point. This yields =“1+2(r)^ (1.20) 

and =(1/r)^. (1.21) 

A number of special crack tip elements have been 
4 5 10 

developed , where displacement method has been used. Also 

7 8 9 

hybrid method has been used to develop singular elements ’ ’ . 

These special crack tip elements contain a singularity in 

the strain field at the crack tip equal to theoretical 
1 2 

singularity . One disadvantage of these special crack tip 

Zl 5 

elements ‘ ’ is that they lack the constant strain and the 
rigid body motion nodes. Therefore they do not pass the patch 
test"'^ and necessary requirements of convergence"*^ 


(1.16) 

(1.17) 

(1 .18) 

(1 .19) 


are not 
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present. R.S. Barsoum ' & R.D. Henshell Both independently 

developed isoparametric cjrack tip elements incorporating 

(1/r)^ singularity satisfying the convergence criteria. Both 

of them dealt vjith 8-noded element and shifting the midside 

nodes tovjards the crack tip at quarter points. R.S. Barsoum 

a general curved element of arbitrary shape for both thick 

and thin shells is proposed for linear fracture analysis of a 

through crack in a shell oi- a plate. The element is derived 

from a degenerate 20'-noded solid isoparametric element using 

reduced integration technique. The (1/r)^ singularity is' 

obtained by the same procedure proposed earlier for t^wo and 

16 1 "^ 

three dimensioiial problems , by placing midside nodes near 

the crack tip at quarterpoints . Several illustrated examples 

ranging from classical solutions to practical problems were 

given to assess the accuracy of solution attainable. 

* 1 8 
Some time later triangular and prismatic elements 

were developed (quadratic and isoparametric). They were 

formed by collapsing one side and placing the mid side node 

at quarter point near the crack tip which show to embody 

(1/r)^ singularity of elastic fracture mechanics and O It ) 

singularity for perfect plasticity. In this work a 8 noded 

isoparametric element is taken and one side of the element. 

is collapsed at the crack tip, see Fig. IC*!) and mapped into a 

square } space. The transformations used are ; 
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8 

X = S Ni ( f s‘i'1 ) 
i=1 


8 

and y = 2 Ni ( -f ) . yj_ 
i=1 


( 1 . 22 ) 


(1.23) 


Ki(f ,n) = (1+ ■?■?!) (1 ("i+nii)' 

- (1 + ff i) !'* * (1 

(1 -fl) ■ (1 -n^) (1+tti) d-''!/) -?l/2 


■where (+1) for corner nodes & (u) for midside 

nodes, x^, are the nodal coordinates for the element. With 
this fonriulation it has been shown that the strains have the 


following singularities ", 


_ 

‘dx 



+ A,| 


(1.25) 


whi-re Aq, b^ & A,^ arc independent of r and are constants 


for (9 = constant). 
Similarly* 



B 


M 

0 

0 

ay 

~^r'’ 

r ' 

av 

3- 

. is 

■5x 

T r 



(1 .26) 


and 


O .27) 
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Similarly, 


D 

iTr ' ^ 




+ D. 


(1 . 28 ) 


■where u & v are the t'wo displacement components 


In the year of 1977 R. Jones and R.J, Callinan^^ presented 
a finite element method for determining stress intensity 
factors in a cracked elastic sheet. Special crack tip elements 
are placed around each crack tip. In the special elements the 
stress and displacements are derived from the exact stress 
function i.e. Airy’s stress function while the continuity, 
of displacements across the special element boundary is 
satisfied in a least square sense. Formulation and numerical 
investigation has been done in case of isotropic elastic field 
and a modified formulation has been proposed in orthotropic 
case but no numerical investigation has been done to verify 
the results. 


^ OkQs ctive s of the pr e sent work 

The present work concerns with the cracks in orthotropic 
plates. Attempt has been made for calculating stress in-tensity 
factors varying poisson’s ratio for cen-trally located crack 
and symmetric side notches on both sides in orthotropic 
rectangular and thin plates under uniform tensile stress 
field. See Fig. iC®"') & 1(f). Another case of a single side 
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notch in rectangular ■ and orthotropic plate has also been 
studied under pure tensile stress and shear stress .and also a 
combination of tensile stress and shear stress. See Fig.1(g). 
Stress intensity factors have been found out for opening mode 
of fracture as -well as sliding mode of fracture i.e. mode I and 
laode II as described in Fig, 1 (b). Use has been made of a 
special crack-tip element "which is so placed around the'crack 
tip and crack lies on the x-axis in the positive direction. 

In the special element the .= tresses and the displacements 
are derived from the exact stress function (Airy’s stress func- 
tion) w/hile the continuity of displacements across the 
boundary of the special element and other general elements 
has been satisfied in a least square sense. Constant strain 
triangular elements have been used as general elements to 
discretize the rest of the structure, -which substantially has 
reduced the domputational cost. The reasons which justify the 
use of the special crack tip element incorporated in this 
work are the following ; 

(1) The stress function used in the formulation of this 

special element is the actual stress function that exists 
in the vicinity of the crack tip which ensures good 
convergence and the error is very much negligible. 

(2)This special element can be of any size or shape which 
gives the liberty to discretize the rest of the structure 
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as per convenience. A fairly coarse mesh can be used 
which in turn significantly increases the computational 
economy and reduces the effort required for data preparation. 

This special element has not been used so far in 
orthotropic elastic field though this has been used in 
previous works and tested to find very attractive results in 
isotropic case. In the present work, the. application of this 
special element has been extended to the orthotropic elastic 
field in three different cases as illustrated in Fig. iC®), 
and 1(g.)} and in all the three cases satisfactory 
convergence of stress intensity factors have been found. In 
chapter 2, equation (2.87) shows that when is equal to zero, 
the evaluation of stress intensity factor is unpossible 
using the general procedure. To overcome the inherent 
limitation of this special- element a scheme has been developed 
to use this special element to find out energy release rate 
to substitute for stress intensity factor whose evaluation 
is not possible as described above ^with the direct method. 

This scheme has clearly been shown in section 3 .2 . The 
effect of applied external stress and poisson’s ratio on 
stress intensity factors and energy release rates have been 
studied in all the three cases. To reduce the computational 
effort and cost advantage has been exploited of symmetry of 
loading and geometry of the plates wherever possible. 
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CHAPTER - 2 
PROBLEM FORM ULATION 


The formulation consists of tuo basically different 
parts. One is the formulation of the special element which 
is placed around the crack tip and the other is the formulation 
of the constant strain triangles which fills up the rest of 
the structure. 


2 . 1 Fgnrul_atlo_n__ o/_ tjae__^^^ strain t r iangles 

2.1.1 Disp lacenient f^ 


■] 

Fig. 2(a) shows a typical triangular element considered 


with nodes isj, m numbered in an anticlockwise order. The 
displacement of a node has two components 3 

^ ^ I . ♦ . 


^i 


( 2 . 1 ) 


( 




The six components of element displacements are listed as vector 



( 2 . 2 ) 


The displacement within an element have to be uniquely defined 
by these six values. The simplest representation is clearly 
given by two linear polynomials 

u = X + cc^y 

V = + “5 X + “57 


(2.3) 
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The six cons "fcCij:! ts ' a > can be solved easily by solving the two 
sets of three simultaneous equations which will arise if the 
nodal coordinates are inserted and the displacements equated 
to the appropriate nodal displacements. Vfriting for example 


(2.4) 



_ a 
_ ^ 

CM 

•^'i 

+ a 

3 

^i 

^d 

- a 

- 1 

, a 

2 

X . 

d 

a 

3 

^d 

^1 

a 

= 1 

, a 

2 

^m 


3 

^m 


We can solve for ^"^5 and in terms of the nodal displacements 
^i’ ^0 ^ obtained finally the horizontal displacement 

u = '2^^ iaj_ + b^ X + Cj^ y) 

u^ + (a^ -i-b^ X + 0^ y) 


d d 




(^.5) 


in which 

a,.= X. y,^ “ X. 


1- "O -m -m ^3 ’ W = ^ ••^2.6) 


and A is the area of the triangle and 


2A= .dot 


1 Yi 

"'^3^3 

Ym 


= 2 (area of the triangle ii^m) 

... (2.7) 
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The equation for vertical displacement are 


V = 


1 


= (a. y) v^+(aj-»bjX+Cjy) 

... ( 2 . 8 ) 


Now we can write in general way 

U =jv|~ = ( I I Njs I ) a® where N is the 

shape function and l is a two by twt- identity matrix. 

Ni = ( X + c^y) /2A 5 

Nj = '-.a^ + bj x+ Cj y)/2.:Vand '^2.9) 


The chosen displacement function automatically guarantees 
continuity of displacements with adjacent elements because the 
displacements vary linearly along any side of the triangle and 
with identical displacement imposed at the nodes, the same 
displacement will clearly exist all along an interface. 

2.1 .2 


The total strain at any point within the element can be 
defined by its three components which contribute to- internal 

'I 

work . 


e =<*'^ y y 

j i 
( 6 xy I 
V / 



('^ X y 1 9 /ax 

0 

9 /9y 

Substituting equation (2.8) vje have 

= B .a'^" = ( B. B . B ) \ 

'' 1 0 m V 

^ { 


u 


u 


L. U 


(2.10) 


V 


••• ( 2 . 11 ) 
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With, a typical matrix given by = , LIN 



aN^/3X 

0 


h 

0 



= 

0 

aNj_/dy 

1 

" ‘2^ 

0 


« • « 

(2.12) 


3Nj_/ay 

aN^/ax 


c . 

1 

^i 





■J 


- 




In general B 

= LN 





(2.13) 


Therefor 

e, the full B matrix look like ; 






0 

aNj/ax 

0 

a\/ax 

0 

B = 

0 

9N^/dy 

0 

aN^ay 

0 




8Nj^/ay 

6N^/ax 

aNj/ay 

aNj/ax 

aNjj^/ay 

aNjj^/ex' 







• • « 

(2.14) 


soj (B )is a (6 x 3) matrix ■'-vhich is actually a strain shape 
function matrix. 


2.1.3 ^ast icity Matrix 

When the elasticity matrix is evaluated in terms of 

20 

principal material direction, it takes the shape given below » 


\ 


Ui 

^2 


^ 1 
/ 

: ^11 

CvJ 

a 

0 ^ 


r ', 

. '"-12 

‘^22 

0 


0 

0 

^66 : 



♦ • f 


(2.15) 



22 . 


•where Q 


12. ^2 


11 - 1"i2i2':?^T ' 


^22 -- " ^^2 


12 . 16 ) 


When we calculate the elasticity matrix -with reference to any 
arbitrary pair of axes (x-y) which makes an angle 'iB' with 
the material principal axes, see Fig. 2;,b), following , is the 
elasticity matrix 

(Q) = (T (Q) (T ... (2.17) 

wh-re ( Q )is 'the modified elastioity matrix and (t ) is the 
transformation matrix. ' '■ 


■ ■ 2 2 

cos b sin © _ 2 sin ©. oos © 

.2 2 

sin © cos G ”2 sin ©. cos © 

I -00 

,-sin©. Gos© sin©. cos© (cos © -sin^" ©) 


( 2 . 18 ) 


Then the stress-strain relation in x-y co-ordinate system is 


f ^) 

(Wty/ ! (Exy^ 




26 ■ ) y(--- ^ 2 .^ 9 ) 


Q'^2 '^22 , ^^26 

^16 °26 ^66 ^ 


where 


cos © + 2 + 2 Qgg) sin © cos © + ©22 sin © 


(Q..|.^+Q22“^ ^ 66 ^ sin 2 © cos2©+q^ 2 (sin © + oos ©) 
Q.^..] sin^ ©-h 2 (q^ 2-J- 2 Qggy sin^ © cos2e + Q^p cos^ © 
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A 6 “ 

“2 

^66^ sin B.cos^©+(Q^2-Q22+2Qgg)sin^©.cose 

~26 ~ 

-2 

^66^ ©.cos© +(Q^2“Q22'^^’'^66^® ©.cos^© 

Q56 ~ (Q']']-rQ22 

-2Q 

2 2 Zi / 

12"^^66'' ©.cos +Qgg(sin^©+cos^©) 


U.20) 

2 . 1.4 

The stiffness matrix of the element ijm is defined from 
the general relationship"' 

^iQ ^/( ) A Q )(Bj_ ) . t dx dy * • • (2.21 } 

where *t' is the thickness of the element. For constant 
strain elements B is independent of x and y i.e. the strain 
is constant throughout the domain of the element^ so, the 
elements of the stiffness miatrix take the following 




where 


is tue area of the tr'i angle 


ijm i.e. A = dx dy 

A 




( 2 . 22 ) 


♦ ♦ » 


( 2 . 23 ) 
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2 . 2 F grmul at i op of sp ec ial c r ac k lemen t 

2.2.1 General_ fj^rmulatlon in _^^otropic_^ 


Consider a polygonal element -with its centre at the 
crack tip see Fig.2(c},. and let the element co-ordinate 


systems he placed with its origin at the centre . Then in polar 

1 y 

co-ordinate system the stress function may be expressed as; 






2n-3 

^2B+T‘' 


.cos 




Here r is the radius of the c ire ascribing circle. See Fig.2(c) 
the number of terms i.e. N, considered in the series expre- 
ssion for cp must be such that it gives accurate values for 
the stresses in the region r < r^ . For any given structure 
as r increases the value of p should be increased. Here, the 

O 
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coefficients d^, ..... ^2’ ...... C2^ are 4XN 

number of degrees of freedom associated with the polygonal 
element. However, the element has an additional three degrees 
of freedom., those are, u^, and w for the rigid body trans- 
lation in vertical and horizontal direction and the rigid ■ 
body relation about the origin. 

The strain energy of this polygonal element in plane stress 
case is given by the following expression 
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Now, any general term of the element stiffness matrix K®iJ 
conceptually means the contribution in the force of the ith 

node due to the displacement at the Jth node."'^ 

Now , 

^^ij ^ jj\ 1 0 


6(p, 




i- //( 2V- ? . <P. .2(1.1) ) (2 (1 ^ ^ j 


1 9^ <P i , 0 


V y 

r ■ 


T& r 5F ae ^ a r^ d 


a 2 ffl . 2 


a 2 <p. q 2 cp 


where , S/ 


^ -TT 2 >) 

0 r^ a r^ 

r dr 

d© 


(2.30; 

a ^ ‘^i 1 

d(p. 

1 

a ^ ^i- 

(2.31) 


”^’r~ ■*' 

d. 

r 



and for i'4 2. N 
i-1 


(0 


1- = v-lj 


r V2-f-1 ^ j _2 




.i/2+1 


^Ti‘2' (i/2-:-1 )9”Cos (i/2->1 )0 

for odd i ..(2.32) 


/A , i/2 

■' j_ = (“1/ — ”"”£7^ (cos ( i/2+1 )©-c os i/2-1 )0) for even i 


(2.33) 


vjhilst for 2 i < 4 N 
i“1 ... 

, .j -rT’ i/2-N+1 

i ^ ' TJYW (i/2-N+l)© -sin(i/2-N-1/©) 


for odd i 


(2.34) 
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i /p„1\j „i/2'“N-!-'1 •; 

and (p '-*1 j == — i72"-^r" (i/2~N+1 )©-sin( i/2-N-l )©) 

for even i . . (2 .35) 


for 4N^i,0 <4 n+ 3, i.e., for rigid body motion is 

ide nt ic ally zer o . 

In practice to find out ■ the integration is done 
numerically, and the shape of the element is largely governed 
by the ease of coupling the remainder of the structure to this 
special element. When the number of sides along the boundary 
of the special element is large then error in approximating 
this polygon by a circle is small. For example in case of 
10,15 and 20 sided polygon the difference in area between 
these polygons and the corresponding circumscribing circle 
is 6.9, 3.0 and 1.6 percent respectively. Consequently 
when the number of sides is greater than say, 20 the numerical 
intigration may be approximated by integrating over the 
circumscribing circle and multiplying the result by the 
ratio of the area of the polygon to the area of the circle. 

27 

Because of the singularity at the crack tip Papaioannou 
22 

and Wilson have found that when only a few terms in the 
series for are retained (e.g. h=1 or 2/, so that is by 
necessity very small, then a very fine mesh is required in 
the vicinity of the special element. This is natural enough 
since in this region the stresses are still changing very 
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rapidly,, Consequently the mesh required for the entire 
structure is quite fine. As observed in it is therefore 
necessary to have r^ quite large , in relation to the crack 
length, so as to escape this region of rapidly changing 
stress. This requires that more terms in the series expansion- 
for be retained and the value of N be greater than that 
previously considered. 


;vith the stress function as given in eqn. (2.24;, the 
expression for the radial and tangential displacements 
U(r,©) and V(r,©) become”^ ^ ; 


N 


U,r,e) = 


1 

YaI 


2 ( “1 ; ) 




(d 


n=1 


2n~1 


2n-3 


cos 


/ 2 n +'1 


}© 


(3.5“h“4cr-} cos )9 C 2 ^_.| sin(^|ii)G + 




ana 


(3.5-n-4cr) sin (^^)© +(-1)*^"''' )" (d 2 ^ (n+1). 

cos (n+1)© +(3-4r3r-n} cos (n-1) © +?2n 


,n+1 /r 


n 


(n+1) © + (3-n“40^, sin (n-1) © 


(2.36) 


N 


V 


(r,©) = -^ 2 (- 1 )^ (^) 


n-'-3 


(c 


2n-1 


n=1 


(2 .5+n"4 Ch')cos 0 ■'‘^2n“1 (2 .5+n-4 ij')sin(“^'^— )0 


cos(^^)© 


.•.,/2n“3' 
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2n-3 „ • ,-2n+1'.cv ^ ^ nII 1 /r .. ,■ , , g^- 

— 2“*“ sin j + (.“1; ■; (d2j^ CnH-3- 

sin (n-1; © ■' (n+1) sin ,n-i-1 j© +'^2n vn-1 ;cos(n4-1 )©-(n+3 
- 4cr) cos (n-1) © ) ,.. (2.37) 

where, cr=;j>for plane strain ando'= l>/(1+2)) for plane stress 

At the ith node, the cartesian displacements are and 

are related to U(r,e) and V(r,©) hy the following formulae ; 

Uf = U^r^,9^) cos ©^ -V(rj_,©^;sin ©i+u^-Wy^ (2.38) 

= U(rj_,©^) sin V(rj_,©j_) cos (2.39) 

where r^. and ©^ are the polar co-ordinates of the ith node 
and are the cartesian co-ordinates of the ith node. 

Substitution of the expression forlj(r,©} andVt^r,©), i.e. 
eq, (2.36) and (2.37) into eqns . (2.38) and (2.39; results 
in a matrix equation of the form L \ = 6, (2.40) 

where, 6 =(u^ , v,| , ^ 2 ^ * '^m number 

of nodes in the special element boundary and L is a trans- 
formation matrix of dimension 2m .x (411+3) . Since it is 
conceivable that the special element may be .coupled to the 
rest of the structure at more points than there are degrees 
of freedom use of least square technique is made to minimise 

discontinuity of the displacement across the boundary of 

i R 

the specie,! element. This yields . 
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\=( L)"”^ L^6 


(2.41) 


so that -when the nodal displacements and are considered 
as the degrees of freedom then the element stiffness matrix 
he comes 


K 


( ( L 


T 


L 


T vT 




( L 


L ; 


-1 


(2.42) 


This formulation of the stiffness matrix is ^/jhen over, it 
can he assembled to the global stiffness matrix. 


2.2.2 Extension of the formulation in orthotropic elastic 

TrelT 


di = 


B 


9y' 




9 2 CP 




xy 


9 ^ cp 
9 X a y 


(2.43) 


where in this case the stress function satisfies the 
compatibility equation ; 


■1 9 ^cp ^1 ^^2 . a ^ cn 1 a ^ fo 

r-- Tx'f " " Tq" 7 " ET aT' ' 


(2,44) 


In an isotropic case where » '^'\2 



and G 



9 


the compatibility equation (2.44) reduces to biharmonic 
equation ; 


^(P=0 (2.45/ 

which which the solution is given by equation (2.1). In the 
general case equation (2.44) can be factorized in the follow- 
ing form 



• • • » 


0^02 0 ^ 0 ^ =0 . . 

•where , for R.rri ,2 , . . .4 

designates the following operation ; 

a 8 

\ ay “^R “8rF' 


31. ' 
(2.46) 


(2,47) 


As shown in 
A ^ 


25 


+ (w 


u 


l_ 

12 


are the roo-ts of the characteristic equation 

E, 


2^* 




12 


2 ■^l 

) M. + ff — ~ = 0 

^2 


(2.48) 


Clearly If J4 - (oc + ip ) is a solution to the characteristic 
equation \,2.48) then so is My- where ( a -iP ) 


And so denoting 2 

2-«j — y ^n d “ X ^ y*oe« (2*49) 

then solution to the compatibility equation (2.44) may be 
written as 

= 2 Re ) + (p^ ( Z 2 ) ^ * (2 *50) ' 

where, and (p^ are the two arbitrary functions of complex 
argumen'ts 2 ^ and Zp respectively. 


Let 

us now define ■ 

two 

systems of polar co-ordinates 


and 

(r-p 9 


0 

0 





^1 ^ 


= ^1 

• 1 

» a « 

( 2 . 51 ) 



is.. 






Tp e 

d- 

= Zp 

• < 


( 2 . 52 ) 

Then r^ 

, . 2 

=UX+ ay) 

+( 

py)^)f ••• 

(2.53) 


^2 

= ((X' 

x2 

’ a y) 

-^( 

p 2 

py) ) ... 

(2.54) 
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s. = a tan ( ) 

' a y 


Sq = a tan (•™— ) 

a y 


(2.55) 


:2.56) 


With 5 5 r25 S 2 defined as above the stress function 

can be expressed in a generalized form ; 


= S (r^ '■ cos ^n+1) + a^^ sin^.m-l ) s^) 

n ' *1 

+ r2 (e^^^^cos(n+1 ) S 2 ■. sin(.nH-1 ; S 2 ) . . (2.57) 

In the limiting case \vhen the elastic modulii tend to those 

of isotropic material, the stress function reduces to the 

24 

same as given by M.L. V/illiams . The stresses which are 
known to have the form : 


Cr^ = 2 Re 


2 ^ *^2 
H'?:” 


( 2 . 58 ) 


(Ty = 2 Re ( 


1 d^ 


dz^ -i- 


(2.59) 


u CQ >1 u p 

= -2Re(/t<3^- 35^) 


(2.60) 


now become the follov/ing 


= S n(n-:-1) (r^^" cos (n-1) s^ -i-a^^ sin(n-1 )s^ ) 


(®4n»1 cos(n"1) S 2 + sinvn-1) S 2 ) ) (2.61) 
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0 ^ = 2 ^ n (n-i- 1 ) ( ((a^-p^) a^^_^-^ 2 ):p a^^) cos(n- 1 ) 

-I- (( a2_p2^ ^4ri“^ ^4n“1 )sin(n“1, (( (cc^-^) 

®4n~l ^ ®4n (ii“'^) S 2 + ((a^-p^) e^^-2apya£^^_^ 

sin (n- 1 ; S 2 } ) ... (2.62) 

= "2 n (n+ 1 ) (r^ icos(n- 1 )s^ (a + P ^ 4 ^ ) +sin(n- 1 ) 

n 

^1 ^4n " ^^4n"1 ^ ^ "^2^"^ ^cos(n-1 ; S2 (a ©4^^^ 

•:-P‘S4n)"- >^®4.n-1^ ^ (2.63) 


If the systGiii of cartesian co-ordinates is placed with its 
origin at the crack tip, see ^'ig. 2 ;d) where for simplicity of 
analysis the geometry of the crack is the reverse of that 

given in Fig. 2i,c;, then the stresses (T' and must vanish 

y 

on y = -I- 0. Noting that for y= +0 we have r^=r 2 =x.... and 
s^ = 0; S 2 = 2 Tt ....... And for y=-0 we have s^= 2 tt ; S 2=05 

r ..- “K' — 

^ _ i 2 — ^ . 

Then the reciuirement that <7^=0 yields 

^ 4 n ~1 = "®4n-1 2 7C^n“1} (2.64) 


and sin 2. % ( n -1 j =0 


(2.65; 


From equation t,2,64) it is clear that n=^, 1, 3/2, 2, ....etc, 

2 a • s., 


'4n 


= cos 2 n (n~ 1 ) i,a 


4n 




^4n-1 . 


o • • 


( 2. 66) 
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Substituting the values of and in terms of a^^_^ 

®-4n givsa in the equations (2.64) and i,2.66), we 
1 Q 

finally ob'cain 



^4n-l (^ 1 *^"^' cos(n+1) s^scos 2 (n-1 ) 

O 

( ■ "p' sin ( n+1 ^ S 2 “ 00 s (n-!- 1 /‘S 2 ))+®’ 4 j^ 

( sin (n+1 / -i-cos 2 7t:(n-1; r 2 ^'^^ sin 2 

U+l) ) . . ... (2.67) 





n(n+1) ^ ^os (n-1) s^- 

2apsin (n-1) s^ )■■' cos 2 7T:.^n-1) 

((^ oc^- P^) “4 a^) cosvn-1) S 2 -( 2( a P ^)“|- ~ 2 xP) 

sin ^n-1) S 2 ))-!-a^^ (sin(n-l) s^ (a^ -P^) 

-!- 2 . cos (n-1/ s^) -■- cos 2 7l(n“1 ) r 2 '^*’"' 

(sin (n-1) S 2 ( - p ^)-!-2ap.cos(n-1 ;s 2 ) )) v.2.68) 




n=i5l 


n(n+1 )i,a^^_^ (cos (n-1) s^ “Cos 2 7 T;(n- 1 ) 

r 2 ^ ^ (cos (n-1) ^2 - sin (.n-1) S 2 )) 


H- a^^ sin(n"1) s^+cos 2 7t(n“1) 

sin (n-1 ) ^2 ) ) ... (.2 .69) 
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~ n-i- 1 ~ n(n+1) ^ ( a cos(n-1 )s^ - p sinin-1 ) 

ii — 2 9 * 9 • • • 

)-cos 2 7ti,n-1 , r 2 


( (X cos(n“'1 j S2+si.n(n“1^ ®2C''“p'^^’’ ■*^) )) 
(r^ ( cos\;n-’1; + sinin=1; 


cos 2 (n-1) ( cos 

\ H'^’l ) S2"«'‘ 

sinvn-1) S2) 




(2.70) 

Let 




p=p^ + i P2 and 

q =q^ ■;■ i q2 . . . . 


( 2 . 72 ) 

where 




, A 4 ^ 

p^ - Re( g-™ - 

'iMo 2 ^2 

4 ^-)= ^“-4- ) 

La 

■■ ■ 

.. (2.73) 

Af2 

Pp =rm( v'-.-. ™ 

-"1 

) ... 

^1 ^1 

• 

( 2 . 74 ) 

= "fte( 'g. ” 

1 V a 

^ E 2 (a^^ 

a^'l 2 

" 

( 2 . 75 ) 


= -Im 


(^_ - 1 )= .ii 32 


£2;#^' 


0-2 “■ \ g 

Then expressions for u and v become 


-&■ "2 2 — ^2.76) 

"^1 E2 ( a + 


2 ^ (n+1) (r^^ (p^ cos n -p2 sin n s. ) » 

n =^5 1 j • • 

^2 ^ (n^^l ) ( c os n S2 9 (p^ “ j — 

sin n S2 (2p^ ^ + P2))) +3-4^^ C^']^(P2 ^ 

p^ sin n s^)+r2’^ cos 271: in-1) iP2 cos n S2 +p^ 
sin n S2) ) ) + - wy ... ( 2 . 77 ) 
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^ (n+1 ) (r^ (q^ cos n s^-q^ sin n 

n=^' s ”1 j • • 

cos 2 71: (n--l) ^ cos n S 2 . (q-| ~ q 2 }~sin n S 2 . 
1 2 q^ -:- q2) ) ) + ^ Cq2 '=0® s^+q^ 

sin n s^ ) -r^ cos 2n (n-1 ) (q 2 cos n S 2 +q^ 

Sin n S 2 ) ) )-i-v^ +wx — (2.78) 


where u^, and. w are the three degrees of freedom for the 
rigid body motions as considered earlier also. 

EjLernep t s t if f ne s s^ D 7 atr iyc 


In this case the strain energy is given by the equation 


19 


V. 


h rr ^ ^ 


hff ( 


^ 2 'T' 2 2 

) dy dx U.79) 


^1 


^12 


This integration is done over the whole area of the special 

element. As before, the expressions for stresses are substituted 

in equation \2.79) and differentiating the total strain energy 

y with respect to each elemental degree of freedom. Vv'e 
s 

can write 

K ®. K = .^-.1 ... > 2 . 80 ; 

ax 

T 

w he re X = ^ , a2 5 a ^ .... a^ n”^ ^ n^ ^o^ ^o^ v2» Si ) 

Here bar is used for denoting parameters in case of ortho- 
tropic elastic field. 
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The 


for 




stiffness matrix may he computed as follows 
i j h N j 



19 ^ 


(2.82) 


while . ~ 0 for i, j •> 4 N. 

-^0 

Here ^^i, (Tyi and "^xyi denote the coefficients of the ith terms 
of the series expansions of che stresses (Tx^'^y and '^xy as 
given hy the equations ’’ 2 . 68)5 ( 2 . 69 ) and ( 2 . 70 ) respectively 
i.e .5 the coefficients of a^^ for i=1 , 4i\i. As discussed in the 
case of isotropic elastic field, it is desirable that the 
special element je connected to the rest of the stmicture 
at more points than there are degrees of freedom. And the 
discontinuity in displacement across the special element 
boundary is minimized in a least square sense. Substituting 
co-ordinates of the nodes in the equations ( 2 . 77 ) and ( 2 . 78 ) 
again results in a matrix equation of the form ; 

(L) = (d} ... (2.83) 

„ T 

where d — ( u^ , v^ , ^ ^ 2 ^ •♦«. ? "'^m ^ C ^ • 84 ) 

is the nodal displacement vector where m is the number of nodal 
points on the boundary of the special element. L is a trans- 
formation matrix of dimension 2mx(4N+3). So the stiffness matrix 
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for the special element, -when the nodal displacements are 
treated as the degrees of freedom, becomes"^ ; 


(K) = ( 


CL)T 


-1 ,'T-^T 


(L) ; (r 


(L) ) 


-1 


(L)" 


(2.85i 


Then (K) is assembled to the global stiffness matrix. Once 
the nodal displacements are solved ~ can be evaluated by 
the relation ; 

\ = (L^ L)”'^ 6 ... (2.86) 


The stress intensity factors and K^j may be evaluated, 

IQ 

once -ciie vector \ is knovm using the formulae " ; 

K^- == LiinO-y (2Tt 

X“^ 0 
y-i*- 0 

= 4 (2 It ... (2.8?) 

v;here a^ is the first element of the vector X. 

Tljjy: Llm ^xy (2 it (x^-i-y-))^ 

X"«* 0 
0 

2 

= 4 2 cxa2) . (2Tt)~ ... (2.88j 

where a 2 is the second elem.ent of the vector 

In calculating Kj and Kjj use is made of the series 
expansions (2.69) and (2.70) for *7^ and ^xy respectively. 
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CHAPTER - 3 
NUMERICA L SCHEME 

3 . 1 Prqgrgijmne_ j-g, P t i on 

The -whole programme is basically an assemblage of twentyone 
sub programme segments meant for different specific operations. 
The main programme is the monitoring part which decides the 
actions and work flow line of the rest twenty programme segments 
or subprogrammes. There are eighteen subroutines and two 


functions as follow 





1 . SICTRX 

8. 

TR0MB1 

15. 

F04AAF 

2, miRX 

9. 

TROMB 

16. 

SIF 

3. ESM 

10. 

SESM 

17. 

STRESS 

4. ASMBLE 

11 . 

MATINT 

18. 

POLAR 

5. SPCL 

12. 

ASMBLS 

19. 

FI 

6. EU4TRX 

13. 

DSMBLE 

20. 

F 

7. UMmx 

14. 

DSFiBLS 



In addition to 

the 

control 

or monitor 

part the main 


programme consists of two more portions which calculate energy 
release rate and stresses in each element after the displacements 
at all the nodes are known if control options for their 
evaluation are given in the data as directed in the beginning 
of the programme listing (Appendix l). A brier account of 
purpose and function of each subprogramme is given below 
in order. 
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1 . 

This subroutine calcxilates the stiffness constants for 
the material to relate stress and strain at any point in the 
domain. It takes material constants like ^2 

as input and after calculating stiffness constants as described 
by the equation(2 .20) it passes these ccnstants to the main 
programme and they are stored. 

2 , 

It calculates (B; as given in equation(2.l4)for each 
general element and passes it to *ESM’ to calculate element 
s ‘tiff ness matrix. 

3. ESM ; 

It calculates the stiffness matrix for each general 
element. It takes vB) of equation(2 .14)and (q) of equationC2 .20) 
as input and calculates each element of each element stiffness 
matrix according to the equation(2.23) . 

4. j 

For each, general element, once the element stiffness matrix 
is formed by *ESK’ it is passed to this subroutine. The function 
of this subroutine is to assemble each element stiffness matrix 
in the global or master stiffness matrix at relevant positions. 

SPCL ; 

This is a control subroutine for evaluation of special 
element stiffness matrix by monitoring the subroutines which are 
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employed to perform this evaluation. Its function is to develop 
the element stiffness matrix for the special element around 
the crack. tip. This calls directly three more subroutines to 
do that. These subroutines are ; 'ELHTRXS 'UI''ITRX* and 'SESM' . 

.6. Eli'IW^j, 

This subroutine develops the L matrix of the equationv2 .83) 
when it is called to do that by SPCL. To do that it takes the 
help of equations (2.77) and (2.78) v^hich describe the total 
displacement field in terms of nodal co-ordinates of the 
special element and the unknown constants i.e. ^?^,^of equation(2.83) • 

7. 

This subroutine generates the elements of stiffness 
matrix (K®) using the equations (2.82). This subroutine 
has to perform a double integration over two dimensional space 
of a complicated function. This is done numerically, To do this 
it calls 'TR0MB1’. 

8 . 

This does the outer integration of the equation (2.79) 

over the whole area of the special element. This subroutine 

integrates IB for the total range of x vdiere IB=/ F(x 5 y; dy . 

The_ limit .for y at any value ofx is a function determined by 

the shape of the element. This subroutine calls F1 to evaluate IB 

at every step of x. The integration is done by ’repeated interval- 

27 

halving and Romberg integration’ technique 
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9. £1 ; 

This function calls another subroutine called ’TROHB' which 
performs the integration of FCxjy.. with respect to y to 
evaluate IB . 

10. TRQivIB 

This does the integration of F{x,y) with respect to y by 

27 

’repeated interval-halving and Romberg integration’ technique 

1 1 . F 

This is the function body which evaluates F(x,yj wioh 
input of X and y values . 

12 . 

This subroutine is developed to perform a number of matrix 
operations like transposition, multiplication and inversion 
to evaluate (K) of the equation (2.85) • Except the inversion 
all other operations are done in this subroutine 'and for inver- 
sion it calls another subroutine called ’HATINV’. It is 
Only here, that final element stiffness matrix for the 
special element is formed. 

13- MA-TEW^ ; 

This subroutine has been developed for matrix inversion by 
28 

partitioning 
14. 3 

Once the element stiffness matrix for the special element is 
formed this is passed to the subroutine ’ASMBLS ’ which assembles 
in in the global stiffness matrix in relevant positions. 
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15. gSMB^^ ; 

This subroutine is needed only when the calculation of 
energy release rate is performed. To he precise what it does 
is just the opposite of ’ASMBLE’. This removes the particular 
contribution of any element stiffness matrix of any general 
element from the global stiffness matrix, when the stiffness 
matrix of an element unaergoes some change, this routine is 
used to remove the previous element stiffness matrix from the 
global stiffness matrix and the new stiffness matrix for che 
element is evaluated by 'ESTI* and assembled to the global stiff- 
ness matrix by 'ASMBLE*. This operation will be discussed later 
in this section only. 

16. DSMBLS ; 

This is a subroutine whose function is identically same as 
'DSMBLE* but this only deals with the element stiffness matrix 
of the special element. 

17. Fj:.4MF ; 

This is a standard subroutine for simultaneous equation 
s olving in the NAG library. Use has been made of this subroutine 
to find out the displacement vector after the force vector and 
the global stiffness matrix have been formed of the equation 

=(K) ... (3.1) 
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Here ^Fj is the force vector, is the- displacement vector 

i.e. the solution vector and (K) is the global stiffness matrix. 

18. ; 

Once the displacement vector is known this subroutine calcu- 
late the stiffness intensity factor by using the equations (2 .86) , 
(2,87) and (2.88) respectively. The equation . (2 .87) is used 
to find out stiffness intensity factor in moae I or the open- 
ing mode and (2.88) is used to calculate that in mode II or 
the sliding mode. 

19. S TRESS ; 

When some points in the special element are given with 
their coordinates in the data, this subroutine calculates 
the stresses at those points if the control option is given for 
that as described in the beginning of the programme listing in 
Appendix I. After the unknown coefficients in the equations (2.67) , 
(2.68) and (2.69) are knovm, this subroutine makes the use 
o f same equations to find out the stresses . 

20. POLAR s 

If the control option for this is given in the data as 
described in the programme listing, it converts all the 
stresses in polar coordinate in the output. 
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S_c for _e ner gy rel eas e rate 

It already has been mentioned in section 1 .4 that ffom 
equation (2.87) it is clear that when is zero, the calculation of 
stiffness intensity factors using equation {2.87} is impossible. 

To overcome this drawback another scheme to determine the energy 
release rate has been introduced in this work. The technique 
that has been used here to evaluate energy release rate is very 
efficient and economic too. The details are presented in 
ref. (29,30). For the sake of completeness the method is briefly 
discussed below ; 

Thi§ method involves calculations based on the direct 
evaluation of changes in the potential energy content as the 
crack progresses. For instance, the potential energy could be 
found for two different positions of the crack tip. In Fig.2(c j 
it has been shown somewhat crudely, a finite element idealiza- 
tion of a structure including a crack. The energy is now 
evaluated for two different positions of the crack tip separated 
bv Aa and approximately we obtain 

Energy release rate = G = ^ — ' (3.2) 

Such approaches were first suggested by Dixon and Pook and 
followed by others . 

The t\'JO separate solutions implicit in such method are 
uneconomic and the direct determination of 'G’ from a single 



solution would be preferable. However, with a simple modifi- 
cation suggested by Parks^^ and Hellen^'^ it is possible to 
avoid such ’double work' . 

To describe this method the system illustrated in Fig.2(c ) 
will again be considered. Let K, a and f correspond to the 
stiffness matrix, displacement vector and the external load 
vector respectively with original position of the crack. Further 
letAa and Ak be the changes in these quantities due to the 
extension of the crack byAa. With original crack positions 
we have 

Ka f = 0 .... (3.3) 

We can now write the change in potential energy due to the crack 
extension as 

i (a+Aa)^ (K+AK) (a+Aa) + la+Aa)^f“ i a\ a -a^f..(3.4) 

Neglecting the second order terms we get 

/\ %- -g a'^ AKa ... (3.5/ 

Therefore, G / ( ..a.t) ••• (3*6 

where ’t' is the thickness of the plate. 

It is evident that to determine G only evaluation of 'a- 
by single solution is required. And to calculate Ail the appropriate 
stiffness changes when the crack propagates by Aa, are to be 
calculated. This is done by recomputing the stiffness matrices and 
substituting the previous stiffness matrices in the global stiffness 
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matrix for those few elements only in the vicinity of the 
crack tip whose geometry gets altered by the crack extension 
including the special element ofcourse . 
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CHAPTER - 4 
RESULT AND DISCUSSION 

To verify the formulation of special element in orthotropic 

case as presented in section 2.2.2 alougwith the whole numerical 

scheme and the computer programme, the results obtained in 

this work are compared with those given in ref. 19. The values 

1 9 

of stress intensity factor have been calciolated in using 

the special element in isotropic case as described in section 

2.2.1. But in the present work the stress intensity factors 

have been calculated using the special element in orthotropic 

case with E^ =£2 andi^v |2 =1^^ = "to bring the isotropicity 

in the orthotropic formulation. The value of Kj that has 

1 9 

been compared with that given in in the case of a thin 
plate of dimension 254 mm x 5 O 8 mm. The uniform tensile stress 
of 689 kPa has been applied in Y-direction while the crack 
is 84.66 mm long. The com.parison of values of Kj has been done 
in two cases. In the first case r^ has been taken to be 25.4 mm 
whereas in the second case r. has been taken to be 16.9 mm, 

V- 

where r is the radius of the special polygonal element as 
c 

shown in Fig. 4. 9 The comparison of the values of Kj obtained 

1 9 

in the present work and in alongwith the theoretical value 
assuming infinitely long plate has been presented below : 
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r* Ref, 19 1 Present work for Theoreticf.1 value 
K-.(kPa.M^) N = 16 1 K-r(kPa, M^) 

Kj (kPa. M^) ^ 


25.4 265.5 266.7 


16.9 263.3 264.1 


273.5 


There are five cases in all for which numerical investi- 
gation have been done. The convergence test has been performed 
in all these five cases. It has been successfully observed 
that the values of stress intensity factors are converging 
very fast with the increase of the number of degrees of freedom 
i.e, the number of elements in the vector Ain equation (2.83). 
When this nuimber is increased from 14 to 16 the average change 
in the values of stress intensity factors are of the order of 
0.069^ which is really good enough for any practical engineering 
purpose . 

This has been clearly mentioned in Chapter-1 (page 3) 

that the stress intensity factors do not depend on the material 
properties of the material concerned. They only depend on the 
geometry of the structure, crack geometry and the applied 
stress conditions. But this is true so far as the isotropic 
elastic field is concerned. In orthotropic elastic field 
it is true that the stress intensity factors do not depend 
on the absolute values of , E2 or even G^2 they depend 
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been performed in last two cases i.e. case-4 and case-5, bad 
there been enough computational facilities . But as there is 
limited computational resources available the study of effect 
of poisson's ratio in case-4 and case-5 has been kept out 
of the scope of the present work. The material properties in 
all five cases are taken as follow ; 

= 4.0 X 10® kPa I E2 = 2.0 x 10® kPa | 0^2 = -5 ^ 

N0W9 all' these five cases will be described and dealt 
with in an increasing order as follow : 

CASE - 1 In this case analysis has been performed for a 
cen'tral crack in a thin plate of dimension 160 mmx200 mrn. 

The length of the crack is 40 mm and r„=15 mm. To save the 
computational tim.e and effort as well, the geometric and loading 
symmetry has been exploited. Only one quadrant i.e. one-fourth 
of the total structure has been considered as shown in Fig. 4. 9. 
It has eventually reduced the problem size to one -fourth of 
the original size. 


The structure has been discretized in a fairly coarse 
mesh in the region far away from the crack-'tip and the special 
element. Gradually the mesh is made finer and finer as it 
approached the special element as shown in Fig. 4.9. There 
are 88 elements in all including 87 general elements and one 
special element around the crack-tip. It has 60 nodes in total. 
In plain stress analysis there are only two degrees of freedom 


•I; 










91883 , 
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per node u and v, so, the global stiffness matrix is of the 
size of 120 X 120 and each of the displacement vector and 
force vector contains 120 elements. A problem of this size 
with N =14 and an integration scheme for the special element 
as described in section 3.1 with sind 

27 

described in 'Repeated interval halving and Romberg integration’ 
takes little over 11 minutes of CPU time in DEC-1090 system, 
where N is the number of elements taken in vector A of equation 
(2.83) . 

For this case convergence test has been performed and 
presented in Table 1 . It can be seen from the table that 
N =14 gives a reasonably good convergence of Ep=0.0659^. Ep 
has been calculated to show the percentage variation of the 
value of Kj for any value of N with respect to that for 
N =16. 

The study of effect of applied stress on stress intensity 
factor Kj with different poisson's ratio values has been done 
for a range of 5000 kPa to 15000 kPa with a constant step 
increment of 1000 kPa in the applied stress. It has been 
presented in Table 9. It is seen that the value of Kj varies 
quite linearly with applied stress as shown in Fig. 4.1 to 
Fig. 4.5 which is of course theoretically expected. To save 
computer time this variation of applied tensile stress in 
Y “direction and calculation of corresponding % values has been 
done in the following manner i 
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Once the element stiffness matrices for all the elements 
including the crack-tip special element is done, they are 
effectively assembled to form the global stiffness matrix. 

When this stage of computation is once reached, a number of 
times the simultaneous equation solving is done with different 
load vectors and after the displacement vector is 
obtained each time, the stress intensity factor is calculated. 
This scheme, instead of having a run afresh each time for 
different applied stress saves lot of computational time and 
cost as well. 

The study of effect of poisson’s ratio on stress 
intensity factor K^. has been studied by varying the value of 
from 0.30 to 0.10 with a constant step reduction of 0.05 . Thus 
5 runs have beer taken fori) =0.30, ^=0.25, 2} =0.20, ^)^^=0. 15 
andX)'=0,10, in this case for each value of a fresh run 
has been taken since when the value of changes all the element 
stiffness matrices undergo some change. The variation of Kj 
and (i.e. the first element of the vector ) with respect 
to the value of D^^s shown in Fig. 4.6 and relevant data is 
tabulated in Table 6. The Table 6 has been prepared for 
applied stress in Y"direction i.e.0“y=5OOO.O kPa. It is evident 
from Fig. 4.6 that decreases very slowly with the increase 
of the value of lJ)and increases with a greater rate 
with respect to value ofI>^but this rate gradually dies 
down as the value ofD increases. 
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CASE "2 Here a thin plate -with symmetric notches on both 
sides under uniform tensile stress in Y^direction is studied. 
The dimension of the plate is 160 mm x 200 mm. The length of 
each notch is 20 mm and r =15 mm. Here also exploiting the 
geometric and loading symmetry only one-fourth of the structure 
has been analysed. Only the first quadrant of the structu]re 
has been discretized for finite element model as shown in 
Fig. 4.10. The rest of the problem description in this case 
is identical to that in case-1 . 

In this case also a reasonably good convergence has been 
observed. This has been shown in Table 2. The percentage error 
in Kj value for H =14 is only 0.073% where Ep is calculated 
taking value for K =16 as standard. In other words this 
is the value of percentage variation in Kj value for N =14 
and N =16. 

The effect of applied stress on stress intensity factor 
Kj with different poisson’s ratio values has been studied 
the same way as in case-l . Fig. 4.1 to Fig. 4.5 and Table 10 
are referred to in this connection. 

The effect of poisson's ratio on stress intensity 
factor Kj has also been studied in an identical manner as 
in case-1 . In this case increases very slowly with the 
increase of the value of poisson’s ratio unlike that in case-l. 
as evident from Fig. 4.7 and Table 7. But Kj increases in the 
same fashion as that in case-1. 
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CASE ~5 In this case a thin plate -with a side notch on one 

side under uniform tensile stress in Y-direction has been 

taken. Tho plate is 80 mm wide and 140 mm long. The crack in 

this case is 15 mm long and p is also same as the crack length 

as shown in Fig. 4.11. Due to the asymmetry in the geometry of 

the structure and crack it is not possible in this case to 

of 

reduce the problem. As a matter ^consequence it is a much 
larger problem compared to those in Case-1 and Case-2. It 
has been so discretized for finite element idealisation that 
it has 104 nodes in all and 161 elements including 1 crack-tip 
special element. Here the special element has 18 nodes compared 
to 10 nodes in Case-1 and Case-2. The global stiffness 
matrix in this case is of the size 208x208 and both the 
displacement and external load vector are of the size 208x1 . 

To execute this problem it takes a little more than 22 minutes 
of CPU time and a huge storage area to store arrays of very 
large size. For this reason to run this problem with a limited 
core area of 50 k here in DEC-1090 is very difficult owing to 
heavy page-fault during execution. 

The convergence observed in this case is also very 
much satisfactory as presented in Table 5. The value of Ep for 
N =14 is only 0.048^c which is the least in all the five cases. 
Perhaps this is due to r^ is equal to the total crack length. 
This is seen in Table 3 that for any practical purpose N =14 
is really good enough. 
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The study of effect of applied stress on Kj has been 
performed in the same way as earlier described in Case-1 
description. This study has been carried out with different 
values of poisson’s ratio as in other cases and variation of 
K-r with respect to applied stress for =0.30 top =0.10 has 
been shown in Fig. 4.1 to Fig. 4,5 and the relevant results 
are tabulated in Table 11 , 

The effect of poisson’s ratio on Kj has also been studied 
in an identical manner as done earlier in first two cases. Here 
the variation of Kj with respect to poisson’s ratio resembles 
that in case-1 as evident from Fig. 4.8 and Table 11. 

CASE -4 In this ®se a thin plate with a 15 mm long side notch 
has been analysed under uniform sinear stress in positive X- 
direction. The plate is 80 mm wide and 140 mm long as in Case-3. 
In this <ase the finite element model is same as that in Case-3, 
(see Fig. 4.11). Only the loads those have been applied to 
achieve uniform shear stress are oriented in positive X- 
direction. In this case the fracture is in pure mode-II i.e. 
the sliding mode. Here i.e. the stress intensity factor 
in sliding mode has been calculated for different applied 
stresses as done earlier in Case-1 , Case-2 and Case-3. The 
effect of poisson’s ratio on Kjj could have been calculated 
as it has been done in all first three cases but the compu- 
tational difficulty due to the shortage of storage in disk 
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and limited core area of 50 k as described earlier in Case-3 
problem description, this has been kept out of the scope of 
the present ■work. 

The variation of K^j -with respect to applied shear stress 
has been studied and sho-wn for only =0.15 in Fig. 4. 2. The 
relevant results of this study has been presented in Table 12. 

It is seen in Fig. 4.2a that Kjj variation is exactly linear 
with respect to applied shear stress for different values of 
poisson's ratio, as it is expected to be. 

CASE -5 This case is an extension of Case-4. The problem 
description in this case is same as that in Case-4 but the 
loading conditions here are different .Here the applied loads 
are 30® inclined to the vertical loads applied in Case -3 
as shown in finite element model in Fig. 4.11. The vertical 
components of the applied loads in this case simulates a uniform 
tensile stress field and on the other hand the horizontal 
components i.e. along the positive X-direction of the applied 
loads simulate a uniform shear stress field. As a result the 
crack in this case is a combination of two modes of fracture 
i.e. opening mode and sliding mode. The stress intensity 

m 

factors in both the modes i.e. Kj & K^j have been found 
out separately. 

In this case the variation of Kj and Kjj have been 
studied with respect to the applied stress cFj normal to 
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a plane ■which is 30 degree inclined to the X-axis (See 
Fig. 4.11). Fig. 4.2a shows this variation which is perfectly 
linear as expected. All the relevant results are presented 
in Table 13. 

The study of effect of poisson’s ratio on stress intensity 
factors has been kept out of the scope of the present work 
owing to the same reason as mentioned in problem, description 
of Case “3 • 

SPECIAL CASE The section 3.2 deals with the calcula- 

tion of energy release rate. This case deals with one such 
case where the evaluation of s'cress intensity factors by the 
direct method is not possible. A thin plate made of 3]>bcP251S 

fibre glass/epoxy has been considered. The required data for its 

20 

material properties has been taken from as follows ; 

= 8 . 0 x 10 ^ psi ; E2=2.7x10^ psi 5 0 ^ 2 ='’ and • 

In this c-^se the characteristic equation (2.48) of Chapter 2 
has’ such roots whose real parts i.e. do not exist. 

The calculation of energy release rate for a plate with a 
central crack of this material has been found out. The dimen- 
sion and crack geometry are same as those in Case-1 . The con- 
vergence of energy release rate G with the increase in nxmiber 
of terms of the vector ^ of equation (2.83) has been observed 
and results are presented in Table 14. 
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CHAPTER - 5 
CONCLUSIONS 

5.1 The results of analysis can be cammed up as follows i 

1 . The special crack tip element can be successfully used 
in orthotropic elastic field and this can be incorporated in 
any standard finite element programme for plane stress problems. 

2. ‘Since this element has no constrains on its size or shape, 
the discretization of the rest of the structure is quite 
flexible . 

3. The convergence of stress intensity factors is reasonably 
good. It is felt that 14 or 16 no. of terms in the vectorAi.e. 
so many number of degrees of freedom excluding u^, and w 
i.e. the rigid body transitional and rotational degrees of 
freedom give very good results. 

4. Energy release rate can also be evaluated if necessary 
using the same special crack tip element with the scheme 
described in section 3.2. 

5. In orthotropic case the stress intensity factors depend 
on some ratios of material properties and poisson's ratio. 

5.2 ^ggestions for extension of the work 

Present work can be extended in number of directions to 
study more about cracks in orthctropic plates. Here only the 



60 f 

effect of poisson’s ratio on computed values of stress 
intensity factors have been studied. The dependence on 

TP IP 

X-t JLJ ^ 

other material property ratios like -rr-^ and 7? " can 

^2 ^12 

be studied to explore the possibility that some range of 
these ratios might exist that may be consibred to be 
relatively safe in a sense that stress intensity factors 
in particular modes for those range are quite low. This is 
of practical importance since the material properties of fibre 
composites can be tailored likewise to get those ratios in 
the preferred ranges. An effective optimization of the 
material properties in this light will certainly lead to more 
economic and safer material design. 

The plates analysed in the present work are all 
specially orthotropic, i.e, the material principal 
directions coincide with the reference axes of the special 
crack tin e.lernent. An extension of the work would be to 
modify for the general orthotropic case. The general elements 
used here can be used in the extended case also. In this 
work, the plates analysed are all laminae, which can be 
extended to all sort of laminates like symmetric, skew- 
symmetric etc . 

A three-dimensional special crack-tip element can also be 
developed with the same concept of using theoretically exact 
stress function and restoring the compatibility of displace-" 
ments across the special element boundaries in a least square 


sense . 
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CASE -1 

CONVERGENCE TABLE FOR Kj 


No. OF TERMS 
TAKEN 

N 

STRESS INTENSITY FACTOR 
F0Ri)=0.30 1 

(KPa M*) 

PERCENTAGE 

ERROR 

S 

6 

2227.88 

12.1CP/4 

8 

2436.73 

3.86P/0 . 

10 

2501 .61 

01 .30/0 

12 

2519.86 

0.58^ 

14 

2532.89 

0.065^ 

16 

2534.56 

0,OCP/o 


TABLE ; 1 


CASE 

CONVERGENCE TABLE FOR 


No. OF TERMS 
TAKEN 

N 

STRESS INTENSITY FACTOR 
FOR)) =0.30 1 

Kj (KPaJvp) 

PERCENTAGE 

ERROR 

6 

2724.75 

13.41?^ 

8 

2962 .28 

5.86^ 

10 

3060.77 

2.7396 

12 

3133.15 

0,4396 

14 

3144.38 

0.073% 

16 

3146,68 

O.OCP% 


TABLE ; 2 



62 


CASE -3 

CONVERGENCE TABLE FOR Kj 


No. OF TERMS 
TAKEN 

N 

STRESS INTENSITY FACTOR 
FOR}) =0.30 1 

Kj (KPa.M®) 

FERCENTAGE 

ERROR 

6 

1886.75 

10.23^0 

8 

2008.44 

4 .44% 

10 

2067.08 

1 .65% 

12 

2091 .46 

0.49% ■ 

14 

21 00.75 

0.04^o 

16 

21 01 .76 

0.0C% 


TABLE : 3 



CAEE-4 



CONVERGENCE TABLE FOR 


No. OF TERl^IS 
TAKEN 

N 

STRESS INTENSITY FACTOR 
Fail's =0.15 i 

'“Kjj (KPacM^) 

PERCENTAGE 

ERROR 

S 

6 

892.74 

13 . 71 % 

8 

986,56 

4.63% 

10 

1022.15 

1 .19% 

12 

1030.74 

0 . 36 % 

14 

1033.91 

0 . 053 % 

16 

1034.46 

o.ocyo 


TABLE : 4 
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CASE -3 

CONVERGENCE TABLE FOR & Kjj 


NO. OF TERl^lS 
TAKEN 

N 

STRESS INTENSITY PERCENTAGE 
FACTCR (MODE-1) ERROR 
FCRVa.=0.15i E 

K^-tKPa.M^) P 

STRESS INTENSITY 
FACTOR (MODE -2) 
FOR =0.15 A 
Kjj(K Pa.M^) 

PERCENTAGE 

ERROR 

6 

616.94 

14.01% 

461 .16 

10.84% 

■ 8 

678.99 

5.36p% 

4 9'+. 83 

4 .3y/o 

10 

701 .81 

2.10% 

506.58 

2 .O &/0 

12 

713.50 

0.55% 

514.75 

0.404 

14 

716.94 

0.071% 

516.91 

0.062% 

16 

717.45 

0.00% 

517.23 

0.0C% 



TABLE : 5 




EFFECT OF 

POISSON’S RATIO IN CASE-1 



POISSON'S ^=Rel_M.~l FIRST COEFFICIENT STRESS INTENSITY 

RATIO (i)^) ^ Tmr/fl OF VECTOR FACTOR (MODE -I ) 

= lm|_X^J Kj(KPaM^) 


= 0.30 

= 0.43642 
= 1 .10624 

1 708. 70 

2534.56 

= 0.25 

= 0.40643 

= 1 .11742 

1733.85 

2377.17 

= 0.20 

= 0.37481 

= 1.12857 

1767.10 

2206.61 

= 0.15 

= 0.33971 

= 1.13964 

1803 .35 

2021 .17 

=0.10 

= 0.30073 

= 1 .15055 

1849.37 

181 7. 51 


TABLE ; 6 


EFFECT OF POISSON »S RATIO IN CASE -2 


POISSON’S 

RATIO 

( ^ 12 ) 

ReZ'^J 

/3 = Im (2 3 

FIRST COEFFICIENT 
OF VECTOR 

^1 

STRESS INTENSITY 
FACTOR (MODE -I) 

Kj (L Pa.M^) 

0.30 

= 0.43642 

= 1.10624 

2121 .37 

3146.68 

0.25 

= 0.40643 

= 1 .11742 

2050.50 

2804.21 

0.20 

= 0.37481 

= 1 .12857 

1981 .71 

2474.60 

0.15 

= 0.33971 

= 1 .13964 

1916.98 

2148.52 

0.10 

= 0.30073 

= 1 .15055 

1856.93 

1824.94 


TABI£ s 7 


OF. .POISSON ^ RATIO c/\SE- 


POISSON’S 

RATIO 

(I)/-.) 


= ReC-^O 

= Im 3 


FIRST COEFFICIENT 
OF VECTOR 


0.30 


0.25 


= 0.45642 
= 1 .10624 

= 0.40643 
= 1 .11742 


1416.99 


1514.80 


= 0.57481 

0.20 1607.99 

=: 1 .12857 


--= 0.33971 

0.15 1 695.28 

= 1 .13964 


= 0.30073 

0,10 1791.29 

= 1 .15055 


i 1 ikdi nhH^.'vr lijM^ 
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STRESS INTENSITY 
FACTOR (MpBE-I^ 
Kj (KPa.kt) 


21 01 .76 


2071 .48 


2007.88 


1900.05 


1760.43 


TABLE s ® 



66 . 


EFFECT OF^. APPLIED_ 

STRESS^ IN_ CASErl 

STRESS APPLIED 
IN Y-DIRECTION 

cry(KPa) 

STRESS INTEN- 
SITY FACTOR 
(MODE “I) FOR 
)-)= 0,30 
'^K.jiKPa.M^) 

STRESS INTEN- 
SITY FACTOR 
(MODE-1) FOR 

KjCKPa.M-S) 

5000.0 

2534.56 

2371 .16 

6000.0 

3041 .50 

2845.41 

7000.0 

3548.42 

3319.68 

8000.0 

4055.31 

3793.91 

9000.0 

45.62.28 

4268.18 

10000.0 

5069.23 

4742,42 

11000.0 

5576.11 

5216.66 

12000.0 

6083.10 

5690,78 

13000.0 

6589.96 

61 65.17 

14000.0 

7096.83 

6639.36 

15000.0 

7603.68 

7115.48 


TABLE : 9 


STRESS INTEN- 
SITY FACTOR 
(MODE“1) FOR 
= 0,20 

Kj(KPa.Mi) 

2206.61 

2648.01 

3089.35 

3530.63 

3971 .95 
4413.32 

4854.64 

5295.96 
5737.28 
6178.61 
6619.93 



EFFECT OF APPLIED STRESS IN CASE -•'I 


STRESS APPLIED 
IN Y-DIFECTION 
cry(KPa) 

5000.0 

6000.0 

7000.0 

8000.0 

9000.0 

10000.0 

11000.0 

12000.0 

13000.0 

IAOOO.O 

15000.0 


STRESS INTEN- 
SITY FACTOR 
(MODE -I) FOR 

'^Kj(KPa.M^) 

2021 .16 

2425 .43 

2829.71 

5233.96 

3638 .1 8 

40/+2 ,43 

4446.65 

4850.85 

5255.15 

5659.35 

6063 .59 


STRESS INTENSITY 
FACTOR (mode -I) 

for: 9 ix = 0.19 

Kj(KPa.M^) 


1617 . 51 
2181 .10 
2544.57 
2908.12 

3271 .61 
3635 .13 

3998.61 
4362.12 
4725.66 

5089.21 

5452.77 


TABLE :9(contd.) 



£F, 

OF APPyED 

STEIESS IN CASE “2 

68. 





c** .-fc. ’ « iuw,vt, ' 1 ■*. V ^ If .1 , ^ 




STRESS APPLIED 
IN Y -DIRECTION 
Cfy(KPa) 

STRESS INTEN" 
SITY FACTOR 
(M0DE~I) FOR 
i>^.^_:^0.30 

Kj (KPa.M^) 

' *1 » '*» -0 v.- -r ' kt i.n- -.V. , 

' S 'TRESS INTEN- 
SITT FACTOR 
(MODE -I) FOR 

Kj(KPa.M^) 

STRESS INTEN- 
SITY FACTOR 
(MODE -I) FOR 
;^^_=0.20 

KjiKPa .M^) 

5000.0 

3146.68 

2804.21 

2474.60 

6000.0 

3776.01 

3365.08 

2969.62 

7000.0 

44.05 .44 

3925.95 

3464 .54 

8000.0 

5034.74 

4486.81 

3959.36 

9000.0 

5664.12 

5047.67 

4454.28 

10000.0 

6293.46 

5608.42 

4949.31 

11000.0 

6922.78 

61 69 .34 

5444 .22 

12000.0 

7552.13 

6730.20 

5939.04 

13000.0 

8181 .39 

7291 .07 

6434.07 

14000.0 

8810.80 

7851 .88 

6928.88 

15000.0 

9440.08 

8412.73 

7423.88 


TABLE 

„.L. 12 
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STRESS APPLIED 
IN Y-DffiECTiON 
0-^(i::Pa) 

m * » » V 1 ( , • .<• »” t ,T( .1 ^ 

STRESS ILTSLSITY 
FACTOR (MODE “I j ' 
F0R}3,^= 0.15 

KT(KPa.M^') 

•IL. 

atft'Trat. ijet-.’g.j. . ■ -i - t- "Te-si 

STEffiSS INTENSITY 
FACTOR (MODE "I ) 
F0R3)^= 0.10 

K^(KPaJ#') 

5000.0 

2148.52 

1824.94 

6000.0 

2578.27 

2189.98 

7000.0 

3007.98 

2555.01 

8000.0 

3437.71 

2920.11 

9000.0 

3867.36 

3284,96 

10000.0 

4297.14 

3649.. 94 

11000.0 

4726.77 

4014.94 

12000.0 

3156.45 

4379.91 

13000.0 

5386.25 

4744.94 

14000.0 

601 5 . 94 

5109.93 

15000.0 

6445 .67 

5474.82 


0 (CONTD . ) 
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E|^^£T_pJ7_A£PL;ffiD_^ 




STRESS APPLIED 
IN Y -DIRECTION 
Cry(KPa) 

STRESS INTEN- 
SITY FACTOR 
(MODE -I) FOR 
i),=0.50 

Kyi^KPa-M"^) 

5000.0 

21 01 .76 

6000.0 

2522.21 

7000.0 

2942.49 

8000.0 

3362.90 

9000.0 

3733.21 

10000.0 

4203.62 

11000.0 

4623.94 

12000.0 

5044.32 

13000.0 

5464.67 

14000.0 

5884.93 

15000.0 

6305 .42 


STRESS INTEN- 
SITY FACTOR 
(MODE -I j FOR 

KjCKPa.M^) 

STRESS INTEN- 
SITY FACTOR 
(MODE -I j FCH 
•)) ^=0.20 

Kj(KPa.M'^) 

2071 .49 

2007.88 

2485.85 

2409.56 

2900.08 

2811 .13 

3314.48 

3212.61 

3728.74 

3614.20 

4142.99 

4015.76 

4557.36 

4417.43 

4971 .67 

4818.96 

5385.87 

5220.58 

5800.27 

5622.06 

6214.56 

6023 .74 


TABLE ; 11 
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EFFECT OF 

APPLED STRESS ON 

IN CASE -4 




"*"^"±" 1 ^ 


STRESS APPLIED 
IN X-DIRECTION 
(j- 3 ^(KPa) 

STRESS INTEN- 
SITY FACTOR 
(M0DE~II) FOR 
i;,.0.15 

Kjj(KPa.M^) 

STRESS APPLED 
IN X -DIRECTI ON 
Q-y. (KPa) 

STRESS INEN- 
SITY FACTOR 
(MODE "II) FOR 
D,^=0-^5 ^ 

Kjj(KPa.M®) 

5000.0 

1034.46 

11000.0 

2275.61 

6000.0 

1241 .30 

12000.0 

2482.84 

7000.0 

1448.14 

13000.0 

2689.66 

8000.0 

1655.24 

14000.0 

2896.38 

9000.0 

1362.13 

15000.0 

3103.55 

10000.0 

2068.94 




( !•*!», » •’ r -wrjjirf 

TABLE ; 

_12 




Bj case -5 


EFFECT OF APPLIED STRESS ON 


Kj & Kjj 


73 « 


'm mm m ,*«•►*.)** vjm. ■ 


SIRESS APPLIED 
'30 DEGREE INCLINED 
TO Y-AXIS 

CTjCKPa) 

STRESS INTENSITY 
FACTOR (mode -I) FOR 
I>,=0.15 

Kj (KPa.M^) 

STRESS INTENSITY 

FACTOR (MODE -II) FOR • 
^^=0.15 

Kjj(KPa.M^) 

5000.0 

717.45 

517.23 

6000.0 

860.96 

620.68 

7000.0 

1004.51 

724.21 

8000.0 

1147.86 

827.63 

9000.0 

1291 .47 

931 .07 

10000.0 

1434.97 

1034.56 

11000.0 

1578.44 

11 37 .'98 

12000.0 

1721 .93 

1241 .38 

13000.0 

1865.37 

1344.85 

14000.0 

2008.91 

1448.29 

15000.0 

2152.35 

1551 .72 

TABLE ; 13 


SPECIAL CASE 


CON VERGEMCE TABLE FOR G 


NO. OB’ TERI4S 
TAICSM 

N 

ENERGY RELEASE 
RATE 

G (N/M) 

PERCENTAGE 

ERROR 




6 

0.9538 

9.81% 

8 

0.9021 

5 . 86 P /0 

10 

0.8803 

1 .35% 

12 

0.8699 

0.16^4 

14 

0.8693 

O.O &/0 

16 

0.8686 

0 . 0 C% 


TABLE 14 
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OOlOn c 
00^00 c 

00300 C 
00400 c 
00500 C 
00600 C 
00700 C 
00800 C 
00900 C 
OlOOO C 
01100 C 
01200 C 
01300 C 
"^01400 C 
01500 C 
01600 C 
01700 
01800 
01900 
02000 
02100 
02200 
02300 
02400 
02500 
02600 
:jj02700 
02800 


APpFNnix 

FTNTTF: FLiFMENT PROGRA” for stress IMTESTTY RACTOR and 
release RRiF RVALUATinM a»’D S-^RESS BMALvgTS USING CRA<: 
SREOIRL ET,ehEUT 

MODEsI, jFIMnS STTFNESS INXENSI'ry FACTOR nsjT.Y 

MODFsaj FlUnS STXFNESS INTFMSI^Y FACTOR AMR DOES STRESS 

MnDF=3:FINDS energy RFlgRASC RATE nWT.Y 

M0DFs!4; finds RNRRGY RELEASE RATE AND DOES STRESS ANAL7 
INODEsl :STRESS nUTpttx IS GIVEN TN (R-THFTA) CO-nRDIHA'I 
IM0 De=2:Strfss nuTpuf ts given tn (y-v) cn-npoiNATE 
ISTRS=1:FTNDS stresses INSIDE SPECIAL ELEnEMT AT GIVER 
ISTRS=!0:0'^HRRWISE 

IDSP=! sDISPLACEMENT OUTPUT IS GIVEN (DISPLACEMENTS AT! 
IDSnsDsOTMEPWISE 

DIMENSION CDRD(2,75) ,NODEC3,100) ,A(150,150) ,R(1505,SNi 

dimension A1 If 150, 150) ,R 1 1 ( 1 5A) ,DELK ( 1 50 , 1 50 3 , DELKR C 1 ^ 

DIMENSION CnRDSC2,75),NSTfl5).CST(2.i5) 

dimension 8f3,63,BT(6,3),SC3,.3),AF(6,6),WKSPCE(150) 

C0MM0n/Al/STM(25,25) I 

CnMM0N/A3/NMAX, JHAX I 

COMMON/ A6/ ALP, BFT, EX, Fy,GXY,R»fUXy 

C0MM0N/A7/N0DES(12,13 

CnMM0N/A8/THCK 

CnH«ON/Ai2/Clf253 

0PEMCUN1T=21 ,DEVICE='DSK',FILE='’DATA11'') 
0PENCfJHTT=22,DEVICE='DSK',FILE='0UT22') i 


105 


02«?00 

03000 

03100 

03200 

03300 

03400 

03500 

03600 

03700 

03800 

03900 

04000 

*^04100 

04200 

04300 

04400 

04500 

04600 

04700 

04800 

04900 

05000 

05100 

05200 

05300 

^5400 

05500 

05600 


QOEMCnHTT=b,D»!:viCE=''I)5K') 

,M0OE,IM0DE,r.5TRS,IDSP 

R*^An (21,*) , WNODR , MEL , ML'.’ ,*■'», Ry , py , GXY ,THTR , RNUXY ,THeK 
REAOC?!,^) , f (CORDCI, J) ,T = l ,2) , J=l,MMOr)E> 

REAnC21,»l , ( (MOOEd, J) ,T=< , 3 ) . J=1 , MEL-O 

00 33 1=1,2 
DO 33 J=t,NNOOE 

33 CnRDCI,J)=COROCT,J)/lOO.O 
0 ^ 34 1=1,2 

on 34 J=:l,fjN0DE 

34 C0RnsCI,J)=0.0 

CALL SMTRyCEX,EY,RMnxy,0xy,THTA.S3 
DO 5 T=1,M 
DO 5 J=1,M 
C TYPE »,T,J 

5 ACI,J3=0.0 

DO 10 IELD=1,meL-1 
lELrlELD 

CALT, BMTrx C 1 Eli , CORD , NODE , R , 8T . AREA ) 

CALL ESMCR,RT,S,AREA,AE,THCK) 

CALL AS’^BLECIEL,AE,M0OE,Al 
C TYPE ♦dEL 

10 CONTIHUE 

READC?!,*) NMAX,JMAX,ALP,RET 
READC21,») TFRN,IERMp1 
MS=Mf^nDE-NNSM+l 

REA9C21,*) (NODESCI,!) ,T = 1 ,MN.5M) 
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Mf 


05700 


Hi::AnC2l ,*’) CC008DSCI,J) ,1 = 1,21 ,J=MS.flMnnE) 

05000 


no 35 1=1,2 

O 50 OO 


DO 35 J=riS,NNnDF 

06000 

35 

CORDSCI, JlsOORDSCT, J)/1DC.0 

06100 


IPLI.srJET.jTETjsTELL 

06500 


CaijT, SPCIjriP'RN,IFRNPl,NNSM,IET,, CORDS) 

06000 


CALL ASM 0tS C NMSM , TED , A , COR DS ) 

06400 


DO 50 1 = 1, D' 

06500 

50 

RCI)=0,0 

06600 


DO 55 I=l,NIjfJ 

06700 


READC21,*> ,DNN,XD,YD 

06800 


RCCDN'’-.1)*2 + 1)=XL 

'^06000 


RCCDNM-1)*2+2)=YL 

07000 

55 

CONTINUE 

07100 


DO 60 NO=:l,NR 

07200 


REAnC2l,f ) ,«I,J 

07300 


I=2*C»'I-l)+U 

07400 


RfDsO.O 

07500 


DO 65 JJ=1,N 

07600 


ACI,JU)=0.0 

07700 


A(JJ,T)=0.0 

07800 

65 

CONTINUE 

07900 


A(I,I)=1 ,0 

08000 

60 

CONTINUE 

08100 


II 

e 

a 

^08200 

44 

Rtl(K)=RCK) 

08300 


DO 46 1=1, N 

08400 


DO 46 J=i,W 
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O8S0O 

08600 

08700 

08800 

08900 

09000 

09100 

09?>00 

09100 

09400 

09500 

09600 

^09700 

09800 

09900 

10000 

10100 

10200 

10300 

10400 

10500 

10600 

10700 

10800 

10900 

^11000 

11100 

11200 


46 AllCI, J)=ACI,J) 

CALTj B’04AAFfA, 150, R, 150, W, 1,8.150, WKSPCE,0) 

TVPF *,TFATT, 

IFCTD.5P.SO,0) GOTO 59 

DO 364 

J=I+1JK=J+1 

364 i«?RITEC22,345) I,RCI) , J.RC.T) .K,R(K5 

345 FnRHATC5X,'RC',T3,') = ',F12.5,5X,^R(',T3.'’» = ',F12,5,5X,' 

l')=',Et2.5) 

59 IFCwnnE-2) 69,69,79 

69 CALL SIFCR.NNSM.STFAC.IFRNPI) 

WRITEf22,2343 STFAC 

234 FORMATCRX, 'STTFMEv5S INTENSITY FACTORS' ,F12, 6) 

TYPE *, STFAC 
IECMODE.E0.2) goto 799 
GOTO 393 

79 REAn€21,*),NETjl 

Dn 22 IFLT4 =’'’ETj1,NEL-1 
lELsJELb 

CALL bmTRXCIEL, CORD, node, R,RT, AREA) 

CALD esncp,i*t,S,AREA,AE,TWCF) 

22 CAUL DSMBT.EflFL.AE.MODE.A) 

lELTjaNEDyTETjsTEtiL 
CALL OSNBLS C NNSN , XEL , A , CORDS ) 
on 169 TJsNS.NNODE 

169 CriRDC1.,TJ)=C0RDCl,IJ)-DFLFX 

DO 21 lELLsNELl.NEL-l 
lELsIELL 



I130n 


CALT, ’='MTRXCIEl.,COPD,finDF!,f?,RT,AREA) 

11400 


CALI, AREA, AE,TWCK) 

11500 


CALL A5MBT.EfIFL,AE,W0nii:,Al 

11600 

21 

CONTIMUE 

11700 


DO 170 IJ=NS,NN0DE 

llfiOO 

170 

CnRnSCl,IJ)=CORDS(l,IJ)-DELEX 

11900 


CALL SPCl,CIFRN,TFRMPl,f|MSW,TEL/CORDS) 

12000 


C ALL DSMBLS r , T EL , A , CO^DS ) 

12100 


on 47 Is:l,H 

12200 


DO 47 Jsl,N 

12300 

47 

DRLKCI, J) = A11(I, J1-A(T,.T) 

12400 


on 49 1=1, M 

^12500 


II 

c 

• 

o 

12600 


on 51 j=i,n 

12700 

51 

S0M=S0R+DELK(T,J)*R11CJ) 

12300 


DFLKRCn=5U« 

12'‘00 

49 

Cn.jTiMiJP 

13000 


STJ’IsO.O 

13100 


on 93 K=1,N 

13200 

53 

Sn*«=.5''M-»-Rl 1 (K')*nELKRCJ^) 

13300 


FKI=6nM/C?.0*THCK4DELEX) 

13400 


WRITEC22,t2t) FKI 

13500 

121 

pnRMATCSX, 'ENERGY RELEASE RATE=',E12. 

13600 


TYPE *, FKI 

13700 


IECM0OK.E0,3) goto 393 

^13800 

799 

DO 701 TSLL=1,NEL-1 

13900 


IEL=IELL 

14000 


CALT, R MTRX ( TEL, CORD, NODE, R, ST, AREA) 
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I4t00 


on 751 T=l,1 

14700 


SnMaO.O 

14100 


DO 801 wi-1,3 

14400 


r)n «ai mi-1,3 

14500 


KK5;2*CM1-1)+N1 

14600 


LsCMOOECMI ,TEI.)-1)*2+M1 

14700 

801 

50M=SnM+BCl,KK)»RfIj5 

14900 


.SM(T)=snM 

14900 

751 

CONTINUE 

15000 


DO 051 K=:l,3 

15100 


SnM=«J,0 

15200 


on 901 J=l,3 

^15300 

901 

SOM=SnM+SfK, J)»SN(J) 

15400 


S5CK,)=snM 

15500 

851 

CnNTINUP 

15600 


I>^CTMnDF’,f!'Q.l) GOTO 852 

15700 

C 

WRIT’S C 22, 607) IEL,SS C 1 ) , ISL ,SSC2) , IRL ,SSC 1 ) 

15800 

C607 

FnRMAT(5X,'v5I(;MXC',T2r') = ',S12.5,5X,'SI<;MY( 

15900 

C 

i'TJioyyc',i?,Fi2.5) 

16000 


GHTO ’101 

16100 

852 

CftLT, POtARCTETi, CORD, NODE, SSI 

16200 

701 

CnflTIMUE 

16300 

393 

IPCTST’RS-1) 394,395,395 

16400 

395 

READ(21 ,») ,NNST 

16500 


on R14 KP=:1,NVST 

4il6600 

814 

RFADC2l,t),N5T(Kp),CSTCl,KP),CST(2,KP) 

16700 


CALL STRESSCNST,CST,NWST) 

16800 

394 

CL0SECU’«IT=21,DEVICE='DSK',PIt,Es''DATA11') 





16900 


Cr.U.SECUOlT=22,DEVTCK='D8K'',P’lT,e=''niIT22'') 

17000 


CTjrj,sEC0"i'^=5,nEVics=!'’nsK''» 

17100 


STOJ^ 

17200 


kmd 

17300 

C 


17400 


SniJROOTTME SM-^RX ( EX » SY , RM’^XY , OXY , TFIT A , S 1 

17500 

c 


17500 


OTMEWSinfj S(3,35,C(3,3) 

17700 


RfJUYX=EYfRNtJXY/EX 

17800 


DMT'?=1 .0-RNOXY=»=RWnyX 

17900 


CCl,n=EX/DMTR 

18000 


CCl ,2)=9Nt5XY»Ey/nMTR 

^18100 


C(2,21=Ey/DMTR 

18200 


C(3,3>sGXY 

18300 


CCl, 31=0.0 

18400 


CC2,3)=0.0 

18500 


on 95 1=1,3 

18600 


DO 95 J=l,3 

18700 

95 

CCJ.I)=CCT,J) 

18800 


Ct=rOSDCTMTR) 

18900 


STsSINOCTHTA) 

19000 


CT4=CT»*4 

19100 


ST4=ST»*4 

19200 


CS22=5QRT(CT4*ST4) 

19300 


C.5i3=CT5(ST»*3) 

^19400 


C.531aST*CC'rY*3) 

19500 


H1=2.0f(CCl,2)+2,0*CC3,3)> 

19600 


H2=CC1 ,1)"C(1,2)-2,0YCC3,3) 



Ill 


19700 

19900 

19900 

20000 

20100 

20200 

20300 

20400 

20500 

20600 

20700 110 

20900 
^20900 
21000 C 
21100 
21200 C 
21300 
21400 
21SOO 
21600 
21700 
21800 
21900 
22000 
22100 
^2200 
22300 
22400 


H3=C(i ,2)-CC2,2) + 2.0>i<CC3,3) 

Hll=Hl=^f’S22 

SC1,1)=C(1,1 )*CT4+H1 l+Cf2,2)*5T4 

Sfl, 2 ) = fCCl,l)+CC 2 , 2 )- 4 , 0 *Cf 3 , 3 ))»CS 22 +C( 1 . , 2 )* CST 4 +CT 4 ) 
SC 2 , 2 )sCCl , 1 )*ST 4 +HU*CC 2 , 214 ‘CT 4 

SCI, 3)=H2’»'C631+H3»CS13 
.$C2, 3)=H2*CS13+H3*CS31 

SC 3 , 31 s:CC(i,i)+CC 2 , 2 )- 2 . 1 'nCl, 2 )- 2 .*CC 3 , 3 )'>ff!S 22 + Cf 3 . 31 *fSTi 
on no 1 = 1,3 
on 110 j=i,i 

SCJ, I)=SCT,J) 

RETURN 

END i 

SUBROUTTNE BMTRXC TETj, CORD, NODE, B, BT. AREA ) 

DIMEMRION CORDC2,751,NqdeC3,1001,BC3,6),B'*’C6,3) 
N1=NOOSC1,IEL) 

N2sN0OEC2,IEIj) I 

N3=NansC3,IE!j) 

X 1 .sC 0 R 0 Cl,Ml) ! 

X 2 =C 0 RDC 1 ,M 2 ) I 

XlaCORDCniO) I 

yi=CnRDC2,N!) I 

y2sCORDC2,N2) 

Y3=C0RDC2,N33 

AREA= 0 . 5 » A 3 S( X 2 *y 3 -X 3 *y 2 -.Xl*y 3 +X 3 *Yt+Xl*Y 2 -X 2 *Yl ) 

AlsX 2 »Y 3 -X 3 *Y 2 



22500 


A2=X3*yi-X11'Y3 

226QO 


A3=5X1*Y2«X25Y1 

22700 


D1=Y2-Y3 

22800 


02=Y3-Y1 

22900 


D3sYi-Y2 

23000 


C1=X3-X2 

23100 


C2=X1-X3 

23200 


C3SX2-X1 

23300 


on 20 I=:l,3 

23400 


DO 20 J=l,6 

23500 

20 

BfI,J1=0.0 

23600 


8Cl,i)=oi 

^23700 


B(1,3'>=D2 

23800 


8(1,53=03 

23900 


8C2,2)=Ci 

24000 


8(2,4)=C2 

24100 


B(2,6)=C3 

24200 


B(3.1)=Cl 

24300 


8(3,2)=D1 

24400 


8C3,33=C2 

24500 


3(3,4)=D2 

24600 


B(3,S)=C3 

24700 


8(3,63=03 

24800 


DO 25 1=1,3 

24900 


DO 25 J=l,6 

^25000 


8(I,J)=BC1, J)/(2,0»4REA3 

25100 

25 

BT(J,T5=B(I,J) 

25200 


RETURN 



113 




25300 


END 

25400 

C 


25500 


SU8R0NTTNE P’SM(B,5T,S, AREA , i^E.THCK) 

25600 

c 


25700 


DIMENSION BC3,6),RTC6,3),.SC3,3),AE(6,6),BTSC6,31 

251100 


DO 100 1=1,6 

25000 


DO 10'' J=1,3 

26000 


snM=o,o 

26100 


on 15 K=l,3 

26200 

15 

5NM=SNM+8TCT,K)7S(K, J1 

26300 


8TSC1, J)=SIJM 

26400 

100 

CONTINUE 

^6500 


DO 40 1=1.6 

26600 


DO 40 J=i,6 

26700 


SUM=0,0 

26800 


DO 30 K=l,3 

26900 

30 

SnM=SnM+BTSCI,K)fB(K, J) 

27000 

40 

AECT,J)=S0M»AREAKTHCK 

27100 


RETURN 

27200 


END 

27300 

C 


27400 


SUBROUTINE ASNBDE(IEL,AE,NODE,A> 

27500 

c 

f « !|i « 4: 1; * If: 4! Ic « ic Ic 1 4= f f f f # 1= ic « f )|c left f 

27600 


DIMErlSIOfi AEC6,6),NODEC3,100) ,AC150.150) 

27700 


DO 45 11=1,3 

^7800 


DO 45 Jl=i,2 

27900 


DO 45 12=1,3 

28000 


DO 45 J2al,2 
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28100 
28200 
28300 
28400 
28?00 
28600 46 

28700 
28900 
28900 C 
29000 
29100 C 
29200 
*29300 
29400 
29500 
29600 
29700 
29800 
29900 
30000 
30100 
30200 
30300 
30400 
30500 
.430600 
30700 
30800 


?4=CTl-l)f2+Jl 


l.= CT2-l)t2+J2 

IGaCMODECTl ,IEL)-l)*2-t-Jt 

JG=CNnOE(T2,IEL)-n»2-fJ2 

AfIG, JG>sA(IG,JG)+AF:(M,TO 

CnfJTIMUE 

RETflR^ 

END 


SrJBROnTTMF PQTjAR C 1E’'j ^ CORD , NODE , SS) 

■iliilitt*tt**ii^***^***^*t**^1^****'*****^* 

DIMENSIOH CORD (2, 75) , NODE C 3 , lOO ) , 5SC 3 ) 


NlsNODECl^IFL) 
N2=NaDEC2#IFD) 
N3=NODE(3,IEL) 
XlsCORDf 1,N1) 
X2=C0RDC1 ,M2) 
X3=C0RDC1 ,N3) 
yi=C0RDC2,Nl ) 
Y2=C0RDC2,M2) 
Y3=CORDC2,N3) 
X=CXi+X2+X3)/3.Q I 
RaS0RTCX4X+Y*Y) t 
S2=SIN(2.0#THETA) I 
STlsCSSCl)+SS(2))/2,0 
SSlaSTl+ST2*C2+ST3»S2 
SS2 sST1-ST2*C2-ST3»S2 
SS3=-ST2=I'S2+ST3*C2 


Y3(yi+V2+y3)/3.0 
THETA=ATAN(Y/X) 
r2=sCOSf2.0*THETA) 
t ST2=CSSCn-8S(2))/2.0 


t ST3*S5(3) 
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30<?00 

31000 1000 

31100 
31200 
31300 C 
31400 
31500 C 
31600 
31700 
31800 
31900 
32000 
•^32100 
32200 
32300 
32400 
32500 
32600 
32700 C 
32800 
32900 C 
33000 
33100 
33200 
33300 
^33400 
33500 
33600 8 


HRITEC22,loOO) IEL,SS1 ,TETj,SS?,TEL,SS3 

FnR'«ATC5X,'SG‘'R<',I3, ') = '.et2.5,5X,'S0MTC',T3, ')=',E12.5,5X, 
t'TAOOTCSlI^'ls'^EIZ.Sl 

RSTHRI ? EM!) 

SOBROOTINE 6PCL C IFRN , TFRNP 1 , NMSW , TET.,, CORDS ) 
i^t********^*********^****'*****'************** 

DIMENSin?} CORDS (2, 75) 

CnRMOM/Al/STM(25,?5) 

COMMOM/A3/NMAX, JMAX 
Cnf1M0fJ/A6/Al.P,8FT,EX,EY.GXY,RMUXY 
COMMON/ A7/M00ES (12,1) 

COMMON/A8/THCK 

CALL UMTRX(IFRN,IFRNP1,NNSM,IEL, CORDS) 

CALI. FLMTRX(I®'RN,IFRNP1,NNsM, TEL, CORDS) 

CALL SESM(IFRNP1,MNSM) 

RETORN 

END 

StfBRnUTTNE TROHBKFTl) 

DIMENSION T1C15,15) 

COMMON/ A2/X1 ,Yl,X2,y2 

COMM0N/A3/NNAX, JMAX 

H1=X1-X2 

DO 8 Kl=l,i5 

DO 8 .11*1,15 

T1CK1,J1)*0.0 
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33700 
33S00 
33900 
34000 
34100 
34200 1 

34300 2 

34400 
34500 
34600 
34700 
34900 3 

*34900 
35000 
35100 
35200 C 
35300 
35400 C, 
35500 
35600 
35700 
35800 
35900 
36000 
36100 
^36200 
36300 
36400 


TlCl,t)=CFifX1)+Fl(X2))*Ht/2.0 
DO 2 Mrri^MHAX 
FR1=H1/C2,0*3<’«) 

IWAX=C2»1:M)-1 
on 1 I=!,TMAX,2 

T1 CN+l,t)=Tl CN+l,l)+FlCFLnAT£T)*FRl+X2) 

TlCW + 1 ,1 )“TtC’«,l)/2.0 + Ht#Tl(N;i,l)/C2.0»*’^) 
nn 3 .t=!2,jmax 
NXMJP2 = '^M6X-J + 2 
F0R>TMls4.0»*(J»l) 

DO 3 Nsl^'JXMJPa 

T1 £N, J) = £FORJVlfflCW+l , J-U-T1 £N, J-l ) VCFORvIMl-l .0) 
FTi=TlC7,JMfiX) 

RETURN 

END 

FUNCTION FlCX) 

C0MM0N/R2/ Xl,yi,X2#Y2 
C0MMaN/A3/NMAX,.THAX 
C0M«0N/A4/Tf 15,15) 

UL=Yl+( (Y2-ri)/CXl-X2))*CXl-X) 

SULsO.O 

X3bX 

CALL TR058CUL,SLL,X3) 

Fl=T£2,.JMaX) 

RETURN 

END 
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36500 

C 


36600 


SUBROnTINE '»'R0MRCfJt,,SljL,X3) 

36700 

c 


36000 


C0MM0W/A3/NWi\X, JMAX 

36900 


CnMW0W/A4/Tf 15,15) 

37000 


H=UL-SLL 

37100 


on to 1=1,15 

37200 


on 10 J=l,15 

37300 

i6 

TCI,J)=0,0 

37400 


T(l,l) = CP*(X3,SLT4)4-FCX3,!!L))*H/2.0 

37500 


on 12 M1=J,NMAX 

37600 


FR=H/2,0**N1 

•^37700 


I«AXJ=2**^1-1 

37800 


DO 11 Il=t,TMAXl,2 

37900 

11 

’r(Nl + l,l)=5T(Nl + l,l)+F(X3,CFt.OATCIl)*FR+RLL)) 

38000 

12 

TfNl + 1 ,l)=T(Nl,!)/2.0+H*’rCNl+1 ,15/2.0*»M1 

38100 


on 13 J1=2,JMAX 

38200 


NXMJP2=NMax-Jt+2 

38300 


FnHJMls4.0t'»(Jl-l) 

38400 


on 13 N1=:1,NXMJP2 

38500 

13 

TCNl,Jn = (FnRJHl*TCNl-H,Jl-l)-TCNl , Jl-1 ) ) / CFORJMl-l’,0 ) 

38600 


RETURN 

38700 


END 

38800 

C 


38900 


FUNCTION FCX,Y) 

^39000 

C 


39100 


C0MM0H/A5/K,L 

39200 


COMMON/ A6/ALP, BET, EX, CY,GXy,RNUXY 
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39300 

39400 

39500 

39600 

39700 

39800 

39900 

40000 

40100 

40300 

40300 

40400 

*40500 

40600 

40700 

40800 

40900 

41000 

41100 

41300 

41300 

41400 

41500 

41600 

41700 

*41800 

41900 

42000 


S(;MXC40),SGMYC40),TXyf40) 

IF'CX,'5’Q. 0,1)1 GOTO 81 
GOTO 82 

St ITCY.f^IO.O.O) YsO.noOi 

82 RlsSOPTC CX + AEjO^Y) *♦?+(: 8rT*Yl»»2l 

R2=.590TCCX-Al,D*y)**2+(B»=:r*Yl**2) 

2t=ATAfKC8F:'rHY)/CX+ALP*Y)l 

Z2=ATA0CC-BET»y)/CX-ALPTY>) 

ITCZl.OT.O.O) 21=3,1415926+21 
I'’CT,2.GT.O.O) Z2=-3. 1415926+2? 

A.SM5S=CAl4i>=(t*2-BF:T**2) 

AMB=C ALP»9ET) 

A08=CALP/BET) 

FAS=!4.0*AT,.P*»2 

REMl=K-(Hi'IXCFL0ATCK)/2,01#2) 

IFCREMl.EO.O.O) GOTO 112 

FT=CX+U/4.0?T=FI*2?V=CPI-1,01 

CPNaCOSOC360»V) 

C2i=COSC21*V) »CZ2=C0S(Z2*V3 
S21=SIMCZ1*V5:SZ2=STS(Z?*V) 

RHl=:S0RTCRl*l'T)/RlyRM2=SaPTCR2**I>/R2 

RT=FI»CFI+1.0) 

XXlaCAS«BS*CZ1-2.0»AMB*SZl)?XX2=C2,0>t!ASMB5*ADB-2.0fAMB)fSZ2 

SnMX(F)=RT=t:CRNiYXXl-CPH*RM2’»'CCASMBS-FASlTCZ2-XX2)l 

SGHYCF)=RT*(CZ1*RM1-CPNTRN2*CCZ2-2,0»AD8=I‘SZ2)) 

XX3 = CAl.P<'CZt-RET»SZl3 

TXYfK)=-RT*CRN14XX3-RW2*CPN*(AljP»rZ2+SZ24f2,0*A!;.P*Ar»B+8FTn) 
GOTO 113 
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42100 112 

42200 

42300 

42400 

42500 

42600 

42700 

42000 

42900 

43000 113 

43100 
43200 
*'43300 
43400 
43500 
43600 
43700 
43000 
43900 
44000 
44100 
44200 
44300 

44400 114 

44500 
^44600 
44700 
44800 


FT = CK/4.'));T=FI*2JV=FT-1.0;CPM=C0.SDf360*Vl 
CZlssCnsCZtl'V) jCZ2sCnS(Z2*V) 

SF.l=STNCZt«V) ?SZ2=SXf>i(Z2 + V) 

RMls5nRT(Rl»*T)/Rt ;RM2=.5aRT(R2**I VR2 

RT=FI*CPI+l.03 

XX4=CSZt»ASWBS+2,0»4MB»CZl) 

Sf5MXCF) = RT*(RNl*XX4+CPfJ*RN2»(SZ2»ASMBS+2.0»ftMB»CZ23) 
5f;MYCK)=RT*(RNlKSZi + CPM*R*^2*SZ2) 

TXY(K)=:-RT*fRNl*CReT*CZl + A DPY.OZl ) «RN2*CPN5 ( BET»CZ2 + ft LP»SZ2 ) ) 
RFM2sL-ClPIXC«'rjOftTCL)/2.0)*2) 

IFCREY2.EO.O,0) GOTO 114 

FT = (L1-1.0)/4.0jT=fi»2;v=FI-1 .0 

CPN=COSD(360*V) ?CZlsCnSCZl*V)jCZ2=ChS(Z2*V) 

SZI*STN(Z1»V5 »SZ2=SlM(Z2fV) 

RN1=S0RTCR1**T)/Rl 

RN2sS0RTCR2*4T)/R2 

RT=FI*CFI+1.01 

XXX4=:f ASMRSl'CZl-2.01'AMB*SZi5;XX5=:(ASMBS*2,0*ADB-2,0*AWB)*SZ2 
SGMXCL}=:RT>t'(RNl»XXX4-CPM*RN25CCASMBS-FAS)KCZ2-XX55) 
5GMYCTi)=RI»(CZl*RNl-CPN*RW2»CCZ2-2.0tADB*.SZ2 51 
XXSsCALPKrzi-RET^SZl) 

TXYCIi5=-RT*CRN1*XX6-RN2*CPN»CALP4CZ2+SZ2»C2.0*:AIiP»ADB+BFT))) 
GOTO 115 

FI=:l./4.0?T=FIl'2?VsFT^l,0;CPM=COSDC360*V> 

CZi*C05CZl*V) »CZ2=C05(Z2»V) 

SZlsSTNCZi^V) jSZ2=SI'«CZ2*V) 

RNisSORTCRl**I)/Rl 

RN2sS0RTCR2»*T)/R2 
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44900 


RI=FI*CFI+1.0) 

45000 


XX7=:CSZl*ASMBy’>+2,0*AMB»CZ1 ) 

45100 


SGMXCL)=RT*{Rfai#XX7+CP!}*Rfl2»CSZ2»flSMBS+2,0»AMB»CZ2)) 

45200 


5GMYCTj5=RI*CR^H*SZ1+C»i^*R'J2*SZ25 

4S300 


TXY(Ll=-RT5!tfRVl^CBE'»’»CZt +» l.P*5Z1)-RN2*CPN»CBET»CZ2 + MP»SZ2) 1 

45400 

115 

XX8=CRNUXY/FX)*fSGMXCK)*SGMYCT..)+SGMXCT*)*SGMY(K)) 

45500 


F3SGHXCK)»5GMX(L)/EX-XX8+SGMYf K)»SGMYCl.VFY+TXYf Kl^TXYCTj) /GXY 

45600 


RETURN 

45700 


END 

45B00 

C 


45900 


SUBROUTTNE UMTRXCTFRNrlERNpi ,NNSM,IEL, CORDS) 

46000 

c 


*^46100 


DIMENSION C0RDSC2,75) 

46200 


CnHMQN/ft2/Xl,Yl,X2,Y2 

46300 


C0MMQN/ft5/K,L 

46400 


C0MMQN/ft7 /NODES C 12,1) 

46500 


C0MM0V/R8/THCK 

46600 


COMM0N/R10/nC40,4O) 

46700 


DO 20 K=1,IFRN 

46800 


on 20 l.= l,IFRN 

46900 


S!1M=:0.0?NNSMMts:NNSM-l 

47000 


DO 5 TN=1,FJNSMH1 

47100 


NIsNODESCIN, J) 

47200 


N23N0DESCTN+1,1) 

47300 


Xl=C0RD5Cl,Nn 

*47400 


Y1=C0RDSC2,N1) 

47500 


X2sC0R0SCl ,N2) 

47600 


Y2=C0RDSC2,W2) 


% 



47700 

47^00 

47000 5 

48000 

48100 20 

48200 

48300 

48400 113 

48500 
48600 
48700 C 
48800 
‘^48000 C 
49000 
49100 
49200 
49300 
49400 
49500 
49S00 
49700 
49900 
49900 
50000 
50100 
^50200 
50300 
50400 


C4LTj TROH^lCFTl) 

SHMsSU'^+FTl 

cnHTi^i'jE 

UCK^IjIsSUM^thcK 

CONT’IMOF 

on 113 7sTFRfl+l,IFRNpi 

DO 113 J=TfRH+l ,IFR9P1 

U(I, J)=0,0 

RETHRM 

SMD 

SOBROOTTNE FLMT9XCIFRW,IF8N91 ,NW5M,XET,,rDRDS) 

DTHEN8I0?! CORnsC2,75) 

CnMMOM/ft6/At>P,8Er,EX,Ey,GXY,RMUXY 

CriMMOM/A7/NnDESCl2,l) 

C0MM0W/A9/ELC25,40) 

Pl = CCCAT4P**2)-CHET**2))-RHUXy)/EX 
P2=2.0*ALP*PET/EX 

QtsAl,o/CEY*( CALP**2) + fBFT*f2)'>)~AT,P*RNUXY/EX 
02=-BET#RMUXY/EX-RET/fEYf f CALP4*2)+fBET»*7)l) 
DO 25 I=UMHSH 
NlaWODEBCI,!) 

XlsCORDSCl,Nl) 

Y1=C0RDSC2,N1) 

RlaSQPTf CXl + AliP4'yi)**2+CBE!’*yi )*»2) 
R2sSQRXCCXl-ALP*yt)**24-CBET»Y1)442) 
ZlsATANCCBET4yi)/(Xl+ALP»Yn) 



122 


50500 
50600 
50700 
50800 
50000 
51000 
51100 
51200 
51300 
51400 
51500 
51600 
•*<51700 
51800 
51000 
52000 
52100 
52200 
52300 
52400 30 

52500 
52600 
52700 
52800 
52900 
^53000 
53100 
53200 


Z2=*^tANCC-BET»Yl)/CXl-ATiP»yi)'> 

IPCZI.LT.0.0) Z1=3.1415926+Z1 
IP’C7.2.GT’,0,0) Z2a-3.14159?6+Z2 
TAB=2,a*AT.p/BF:T 
Ks2*I-lsl4=2»I 
ELCK,TFRi'|Plls1 ,0 

EI,CT4,TFRI^o1)s:0.0 I 

nn 30 J1=1,XFBN,2 

FNl=(Jl+1.0)/4.0jIFMle2.0*FNl 

R91=SORTCRi**TFMi3.RriP-.sgRTCR2»*IFNl ) | 

CON=Cn50C360i:CFMl-1.0)) | 

CZ1=C0SCFMI+Z1) ?CZ2=Cns(FMl*Z2) 

SZl=STN(FMi*zi)?sZ2=STNCFMl*Z?) 

FNllaFfn + 1.0?AB=0.0?CR=0.0 | 

XX9=CP1#TAB+P2) . 1 

ARaFNI 15C9'^1»(P1»CZ1-P2*SZ13-RN21'CPM*(CZ2*(P1-TAB*P2)-SZ2*XX9)) 
X10aCQltCZl-Q24SZ!)?Xlla(01-TAB»02) I 

CB=FH11*CRN14X1O+rn2#CPN»CCZ2»X11-SZ2*(01»TAB+02))) | 

ELCK,Jl)aABjELCL,Jl)aCB | 

CnNTIWUE 

DO 40 J2a2,tFRN,2 I 

FV2aJ2/4.0?IF^2a2*FM2 I 

RMlaSORT(Rl**XFN2')?Pri2=5QRT(R2»*IFN2) | 

CZlaCOSCFJ^2*Zt) fCZ2=CnSCFN2#Z?) 

SZlaSIN(FN2*Zl)?SZ2=5l»4CF*J2*Z?) I 

FM2iaFN2+1.0?CPNsC0SDC360»CFN?-l,0)) | 

BRaFN2l*C9Nl*CP2*CZl+Pl*S74 5+PN2*CP>?>t'fP2*CZ2+Pl*SZ2)) | 

DB=FN2i*CRNl*(02*CZl+Qi*SZl)-RN2»CPM*(02*CZ2+01*SZ2) ) | 
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53300 

53400 

53500 

53600 

53700 

53800 

53O00 

54000 

54100 

54200 

54300 

54400 

%4500 

54600 

54700 

54800 

54900 

55000 

55100 

55200 

55300 

55400 

55500 

55600 

55700 

,55800 

55900 

56000 


El.CK,J2)=RBjELCTi,J2)=08 
40 CONTINU!? 

25 CONTINUE 

RETfFRV 
EWD 

C iif^:iitt**lf*‘^*t**1f**'¥***'¥**ilt'¥*^* 

SOBRQOTTNE SESM f IFRNpl , NNSM ) 

COMMO»*/Al/S'^MC25,25) 

COMMO’^/ft9/ELC25,40) 

CnMHO'4/TNV/ELTELXf40,401 

cnMMov/4io/nc40,40) 

CnM'*0V/All/EEC40,25) 

DTMENSin?i ELTC40,25),ELTELC40,40),ET,TEl)TC40,40) 
OIMENSIOM ELTELS(40,40) 

OTMENSIO’’! EET(25,403,EeTU(25,40),TESTC40,401 
DOUBLE ORECTSTON ELTELX 
NNSM2=2*N'?SM 
DO 50 I=1,NNSM2 
DO 50 JslrlERNpi 
50 ELTCJ,I1=ELCI, J) 

DO 60 K=1,IFRNP1 
DO 65 L*1,I’='RNP1 
SUHaO.O 

DO 70 J=1,NNSM2 

70 SUMsSnM+ELT(K,J)4ELCJ,L> 

ELTELCK,L)=SUM 
65 CONTINUE 



1?4 


t 


56100 

60 

contiwuf: 

se'Joo 

C 

SCftTil'^G IS nONE WITH I00f)0.0 

56300 


on 23 IT=1 ,TFRNP1 

56400 


DO 23 JiI=:t,TFRNPl 

56500 

23 

ELTFLSCII,JJ)=EtTFLClT,JJ)*1.0 

56600 


on 111 IS=l,lPRNpt 

56700 


DO 111 .TSsl,IFRMpi 

56800 

111 

El.l’FLXCIS,JS)=EtTFLSCIS, JS) 

56900 


C4LD MATINVCIP’RNP1,T01 

57000 


DO 200 11 = 1 ,IPRNpi. 

57100 


DO 210 Jlsl^IPRUPl 

57700 


SnH=0,() 

^7300 


00 220 I2=i,IFRNpi 

57400 

220 

SUM=SnM-*-EI;TELXCri,I2)*EDTELSCT2,Jt ) 

57500 


TEST(T1,J1 )=SOH 

57600 

210 

CONTINUE 

57700 

200 

CnNTI*?UE 

57800 

C 

SCALING IS DONE HERE 

57900 


DO 230 IH=1,IERNP1 

58000 


DO 230 JHal,IPRNPl 

58100 

230 

ELTELTCTH,JH)=ELTELXCIH,JH)»1.0 

58200 


DO 75 K1=1,IFRNP1 

58300 


DO 80 L1»1,NN5M2 

58400 


SUMsO.O 

58500 


DO 85 J1=1,IFRNP1 

^58600 

85 

SUMaSnM4-ELTELTCKl,Jl)*ELTCJl ,M) 

58700 


EECKl,Ll}=SnM 

58800 

80 

CONTINUE 



58900 

75 

Continuf: 

59000 


on 90 Ilsl,TFRNPl 

59100 


DO 90 J1=1,0NSM2 

59500 


E*^Tf J1,T1)=EECI1 ,J1) 

59100 

90 

contihue 

59400 


DO 95 K2=1,NNSM2 

59500 


DO 100 L2:=1,IPRMP1 

59600 


s?m=o.o 

59700 


DO 105 .12=1,IERNP1 

59900 

105 

SnMsSOM+SETCK2,J2’>*U€J2,L2) 

59900 


EETU(*f2,L2)=SUH 

60000 

100 

CONTIMUE 

^60100 

95 

CONTINUE 

60200 


DO 110 K3=1,NNSM2 

60100 


DO 115 l,3sl,NNSM2 

60400 


SOMaO.O 

60500 


DO 120 J3=1,IFRNP1 

60600 

120 

SnMsSOM+EETUCK3,j’3)*EECJ3,L3) 

60700 


STMfK3,D3)aSUM 

60800 

115 

CONTINUE 

60900 

110 

CONTINUE 

61000 


RETURN 

61100 


END 

61200 

C 


61300 


SUBROUTINE ASNBDSCNNSN, TED, A ,CORDS) 

^1400 

c 

•$1f:)l/iilt^***t*****************^********* 

61500 


dimension 4(150,150) 

61600 


DIMENSION CORDSC2,75) 



61700 


C0MM0M/A1/STMC25,251 

61900 


COMMON/ A7 /nodes C 12, 1 ) 

61000 


00 45 11=1, NNSM 

62000 


DO 45 Jt=l,2 

62100 


DO 45 12=1. ,NNSM 

62200 


DO 45 J2=i,2 

62300 


M=CXl-l)*2+JUL=Ct2~15*2+J2 

62400 


in=CCN0DESCTl,l)-1)*2+Jl) 

62600 


JOaCCNnnES(T2,l)-l)*2+J2) 

62600 


ACIC,.TGl = «i(XG, JG) + STMCM,L1 

62700 

45 

CONTINUE 

62800 


RETURN 

11^62900 


END 

63000 

C 


63100 


subroutine DSMBLS CNNSM , TETj , A , rORDS) 

63200 

C 

H;ili*tt*Mtt*4'$*'¥****t‘***************** 

63300 


DTMENSTON AC150,150) 

63400 


DIMENSION C0RDSC2,75) 

63500 


COMMON/ A1 /STM (25, 25) 

63600 


COMMON/ A7 /NODES (12,1) 

63700 


DO 313 TZ=1,75 

63800 


DO 313 JZ=1,75 

63900 

313 

C0Ri>SfIZ,>IZ)s0,0 

64000 


00 45 Il=t,NNSM 

64100 


DO 45 Jl=l,2 

*64200 


DO 45 12=1, NNSH 

64300 


DO 45 J2=l,2 

64400 


M=(Ii-l)»2+Jl»L®€I2-l)»2+J2 



64500 


IGsCCM0nESCTl,l)-l)=<‘2+Jl) 

64600 


JGsCC!^OOESCT2,l)-t)*2+J2) 

64700 


ACIG,JG> = -^CIG,JG)-STMCM,D) 

64800 

45 

CONTIWIJE 

64900 


RETOR*^ 

65000 


E^D 

65100 

C 

i^Hi:^:$1ftt*V*^********************* 

65200 


SUBROUTINE DSMBUEf lEL, AK,N0nB.Al 

65.300 

c 


65400 


OTMENSION AEC6,6) ,MOOE(3, 1 OO) ,A( ISO, 1501 

65500 


on 45 11=1,3 

65600 


on 45 Jl=l,2 

^65700 


on 45 12=1,3 

65800 


DO 45 J?=l,2 

65900 


Ms(Tl-l)»2+JUU='(I2-l)»2+J2 

66000 


IG=( (NOPEf 11 ,IEL)-1)*2+.J1) 

66100 


JG=((N0nECl2,TETj)-l)*2+J2) 

66200 


ACIG,JG)=ACIG,JG)-AECW,L) 

66300 

45 

CONTINUE 

66400 


RETURN 

66500 


END 

66600 

66700 

C 

f ♦ j*! f f ♦ :|t *#♦♦♦♦ ♦ 

66800 


SUBROUTINE WAVINVCH,IP) 

66900 

C 


^67000 


COMMON/INV/AC40,40) 

67100 


DIMENSION XC40),YC40) 

67200 


DOUBLE PRECISION A,X,y 
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67300 


in = 1 

67400 


IFC M - 1 ) 17,17,2 

67500 

2 

A(l,l)=1.0/A(l ,1) 

67600 


DO 12 K = 2 , M 

67700 


N = K - 1 

67SOO 


DO 11 I = 1 , N 

67900 


yci) s 0,0 

68000 


XCI) a 0.0 

68100 


DO 10 J = 1 , N 

68200 


y(I) S YCT) + ACI.Jl * ACJ.ff) 

68300 


XCI) = X(T) + ACK.J) ♦ ACJ,T) 

68400 

10 

CONTINUE 

*•68500 


ACK,K) = ACK,K) - ACK.I) f Y(T) 

68600 

11 

CONTINUE 

68700 


IECACK,ff)) 14 , 15 , 14 

68800 

14 

ACK.K) = t.O/A(K,K) 

68900 


DO 12 I = 1 , N 

69000 


ACK.I) * -XCI) * A(K,K) 

69100 


ACI.K) = -YCI) * ACK,K) 

69200 


DO 12 J » 1 , N 

69300 


ACJ,I1 = ACJ,T) - YCJ) * ACK,I) 

69400 

12 

CONTINUE 

69500 


RETURN 

69600 

15 

ID = 2 

69700 


TYPE 40, ID 

j^69800 

40 

for«at(5x,'watrix is singular as id in tnv s', id 

69900 

17 

RETURN 

70000 


END 
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70100 

C 


70200 


SnaROUTTNP' SIRCR,NNSMrSTFftC,IFRNpi ) 

70300 

c 


70400 


DIMENSION R(1503, 0(251, Clf25) 

70500 


CnM«0N/A7/N0DESCi?,l 3 

70600 


C0M«0W/A6/AT.P,BET,EX,EY,GXy,RNUXy 

70700 


CnMMOW/All/EEC4n,253 

70800 


M=2*N1VSM 

70900 


DO 555 T0al,NNSM 

71000 


NlaNODESCIQ, 13 

71100 


K=2*N1-1 

71200 


L32*N1 

4?71300 


101=2*10-1 

71400 


I02=2*10 

71500 


CCI013=8CK3 

71600 


C(I023sRCl.5 

71700 

5'55 

CONTINUE 

71800 


DO 111 II=1,IFRNP1 

71900 


SUMaO.O 

72000 


DO 222 JJ21,M 

72100 

222 

SUHaSUM+EECII,JJ3*CCJJ3 

72200 


C1C1I3=SU« 

72300 

ill 

CONTINUE 

72400 


DO 126 Dpai,M 

72500 

126 

WRITE (22, 456 3 DR,CCIjP 3 ,LP,Cl Cl»P3 

^72600 

456 

FORMAT(5X,'C(',T2,')a',F12.4,5X,'Cl(%I2.'3*'.F12.43 

72700 

C 

type ♦, CtCt3,Cl(23,Cl(33,CU43»Cl(53,ClC63#Cl(73 

72800 


SIFAC=C3.0/4,03*SORTC2.0*3.14159263*2,0*(ALP/BET3*C1C13 



72900 
'73000 
73f00 C 
73200 
73300 C 
73400 
73500 
73600 
73700 
73800 
73900 
74000 
^74100 
74200 
74300 
74400 
74500 
74600 
74700 
74800 
74900 
75000 
75100 
75200 
7S300 
^75400 
75500 
75600 


no 


RKTURfJ 

EWD 

S08R0UTTNE STRESS(«Sr,CST,N»JST,TFRNl 

dimension NSTCNNST) ,CST(2,NNS7) ,S0SPXC15) 
DIMENSIOM S<1SPYC15) ,TSPXyC15) 
CnMMON/Ai2/CK25) 

COMMON/ A6/AT,P, BET, ex, Ey,GXy,RMUXY 
ASMBS=CAi:jP>t<*2-BET**2) 
AMBsA!jP*BET?Ar>B=ALP/8ET 
FASa4,0*CALP»*2) 

DO 815 Is: 1,MNST 
SX=0 . 0 ? S Y=0 . 0 ; TY=o . 0 
X=CSTCl,n?Y=C5TC2,T) 

H1=R0RTC CX+ALP*V)**2+fBET*Y)*»25 
R2=v5QRTC CX-ALI»*Y)»42 + CBFT*Y)**2') 

ZlsATANI (BE'’’*Y}/(X+ALP*Y) 1 
Z2=ATAMC-C0ET»Y)/fX-ALp*y)) 

IF-CZl.LT.O.O) Z1 = 3.1415926+Z1 
IFCZ2,CT,0,0) Z2=-3,1415026+Z2 
DO 816 IF=1,IFRN 
IFS*IE"(IRIXC8’LPATCK3/2.0)»2) 
IFCTFS.RQ.O) COTO 817 
FIsCIFtn/4.0?II=FI>l‘2 
yaCPI-l .0) ?CPN3C0SDC360*V> 

CZl*COSCZl*V) :CZ2=COSCZ2>l'V) 

SZlaSTNCZl*V) ?SZ2=SINCZ2*V) 



m 


75700 

75800 

75900 

76000 

76100 

76200 

76300 

76400 

76500 

76600 817 

76700 
76800 
f76900 
77000 
77100 
77200 
77300 
77400 
77500 
77600 
77700 
77800 
77900 

78000 818 

78100 
^78200 
78300 816 

78400 


RHls!SORTCRl*fII)/RljRM2=SORTCP2»»TIl/R2 

RT:s:ri*CFI + 1.01 

XXlaC A5M8S*CZ1 -2.0 )(cAMB1£SZI ) 

XX2=C2.0=|!ASM8SilcAD»-2,0*AMR)*SZ2 

SnMl=ClCIF)*RT*CRHl*XXl-CRN*R»12’<'CCASHBS-FAS)*0Z2-XX2>l 
SnM2=:Ci CiF)*Rl!leCCZl*RMl-CPM*R«2*CCZ?-2.0*ADR*SZ?) ■) 
XX3aCALP«CZl-RET*SZ1 ) 

SnM3=ClCIF)*RT»CRWl*XX3-CPNl'R>r2*CALP*CZ2+SZ2*C2.0*ATiP<!AOB+BF;T3 
(SOTO 818 

FT=CIF/4,0) ?II=2*FI 
VsCFI-i.O) ?CPMs:CQSDC360»V) 

CZ1=C05CZ1*V)?CZ2=C0SCZ21'V) 

SZlsSTNCZl^^'V) ;SZ2aSINf Z2»V) 

RMiaSORTCRl*l'II)/Rl 
RN2aSQRTCR2**II)/R2 
RI=FI*CFI+1.01 
XX4=CSZl=j!ASM8S+2.0KAMB»CZl ) 

XX5aCSZl»ASWBS+2,0»AMB*f!Zl) 

S0M1=C1CIF)*RT*CRN1*XX4+CPN»RW2*XX5> 

SnM2=ClCIF)*RI»CRNl*SZl+CON*8W2*SZ2) 

XX6a(RET*CZl+ALP»SZ! ) 

XX7=C8ET»CZ2+ALP*SZ2) 

SnM3=ClCIF)=»RT»CRNl»XX6«CPM*RM2*XX7) 

SXsSX+SOMl 

SY3sy+sim2 

TY3TY+S053 

continue 

SOSPXCDsSX 
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78500 

78600 

78700 

78800 

78900 

79000 

79100 

79200 


Sf^soycDssy 

T5PXYCI)=TY 

WRITEC22,212i) T ,SGSPXCI) ,I,SOSPY(mT,TSPXYCX) 

2121 FnRMATC5X,'SGSPX(',T2,''> = %E17.5,5X.'SGSPyC',T2,'5 = %P:i2.5 

1,5X,'TSPXYC',I2,')s:%E12.5) 

815 CONTINUE 

RETURN 
END 



